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ABSTRACT 


A  semi-empirical  computer  model  of  the  lower  thermosphere  has 
been  developed  that  provides  a  description  of  the  composition  and  dynamics 
of  the  thermosphere  (Killeen  et  al.,  1992).  Input  variables  needed  to  run  the 
VSH  model  include  time,  space  and  geophysical  conditions. 

One  of  the  output  variables  the  model  provides,  neutral  density,  is  of 
particular  interest  to  the  U.S.  Air  Force.  Neutral  densities  vary  both  as  a 
result  of  change  in  solar  flux  (eg.  the  solar  cycle)  and  as  a  result  of  changes 
in  the  magnetosphere  (eg.  large  changes  occur  in  neutral  density  during 
geomagnetic  storms).  Satellites  in  earth  orbit  experience  aerodynamic 
drag  due  to  the  atmospheric  density  of  the  thermosphere.  Variability  in  the 
neutral  density  described  above  affects  the  drag  a  satellite  experiences  and 
as  a  result  can  change  the  orbital  characteristics  of  the  satellite.  These 
changes  make  it  difficult  to  track  the  satellite's  position. 

Therefore,  it  is  particularly  important  to  insure  that  the  accuracy  of 
the  model's  neutral  density  is  optimized  for  all  input  parameters.  To 
accomplish  this,  a  validation  program  was  developed  to  evaluate  the 
strengths  and  weaknesses  of  the  model's  density  output  by  comparing  it  to 
SETA-2  (satellite  electrostatic  accelerometer)  total  mass  density 
measurements. 

The  SETA-2  instrument  was  flown  on  a  low  altitude  satellite  in  1982 
during  solar  maximum  of  solar  cycle  21 .  The  data  set  used  includes  over 
186,000  data  points  gathered  between  *  7  May  and  22  July  1982.  As  a  result, 
a  large  statistical  distribution  was  available  to  determine  deviations  of  the 
VSH  model  density  from  the  SETA-2  data. 


Density  output  from  the  MSIS-86  (Mass  Spectrometer-Incoherent 
Scatter)  empirical  model  was  also  compared  with  the  SETA-2  data  set.  In 
this  way  comparisons  of  VSH  and  MSIS  deviations  from  SETA  2  could  be 
used  to  test  their  performance. 

A  one  to  one  comparison  was  made  between  the  VSH  and  MSIS 
models  and  the  satellite  data.  The  average  relative  error,  or  percent 
difference,  for  the  models  was  calculated  for  each  point  and  was  plotted 
against  their  frequency  of  occurrence. 

A  significant  shortcoming  in  the  model  was  discovered  as  a  result  of 
conducting  this  validation  study.  By  binning  the  data  in  a  variety  of 
parameters  we  were  able  to  evaluate  latitudinal  and  geomagnetic 
dependences  of  the  two  models.  Initially,  the  plots  showed  there  were 
conditions  when  the  VSH  model  overestimated  neutral  densities.  The  data 
was  binned  by  FI 0.7  values  and  was  discovered  that  the  area  of 
overestimation  occurred  during  times  of  high  FI 0.7.  This  problem  was 
corrected  and  subsequent  results  showed  lower  standard  deviation  values 
by  10%  (for  this  particular  case).  The  study  also  showed  that  the  models 
provided  the  best  standard  deviations  (6-8%)  in  the  low  latitude  regions  and 
during  low  magnetic  activity  and  conversely,  poorer  standard  deviations 
(15-18%)  in  the  high  latitudes  and  during  high  magnetic  activity. 

Atmospheric  density  is  not  the  only  cause  of  drag  that  a  satellite 
experiences.  In-track  winds  that  are  anti-parallel  to  the  direction  of  the 
satellite  can  also  induce  drag,  affecting  the  orbit  of  the  satellite.  In 
particular,  winds  in  the  polar  cap  region  can  reach  speeds  over  1000  m/s,  or 
roughly  10-15%  of  the  satellite  speed. 

In  an  attempt  to  quantify  some  of  these  questions,  neutral  winds 
from  the  VSH  model  were  calculated  throughout  a  typical  orbit.  Many 


different  orbital  tracks  were  used  where  different  parameters  were 
changed,  such  as  varying  the  local  time  "slice"  of  the  satellite’s  orbit  in 
order  to  find  in-track  winds  in  the  polar  region  where  these  circulation 
cells  exist. 

These  winds  are  quite  important,  with  results  showing  that  for  a  742 
m/s  wind  the  satellite  drag  changes  by  23%.  This  is  for  a  satellite  crossing 
the  polar  cap  region  across  the  14-02  local  time  meridian  during  high 
magnetic  activity.  However,  this  is  a  very  specific  case  during  a  certain 
day,  the  circulation  cell  is  highly  variable  in  it’s  position.  Nevertheless,  it  is 
important  to  know  what  local  times  the  convergence  area  of  the  circulation 
cells  occur  and  the  direction  of  these  winds. 
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CHAPTER  1 


INTRODUCTION 


1.1  General 


The  objective  of  this  thesis  is  the  validation  of  the  neutral  density 
output  of  the  VSH  (Vector  Spherical  Harmonic)  model  and  the 
determination  of  the  atmospheric  drag  effects  on  artificial  satellites  duo  to 
in-track  winds  using  VSH  model  neutral  wind  output. 

When  planning  a  satellite  mission  it  is  important  to  consider  the 
atmospheric  drag  effects  on  the  satellite.  Drag  affects  the  lifetime 
prediction,  design  and  control,  tracking,  on-board  fuel  requirements,  and 
the  reentry  (Marcos  et  a!. ,1989).  Therefore,  accurate  specification  of 
neutral  density  is  quite  important  and  one  of  the  primary  objectives  of  the 
VSH  thermospheric  model. 

To  validate  the  model,  neutral  density  measurements  from  the 
satellite  electrostatic  tri axial  accelerometer  (SETA-2)  instrument  are  used. 
The  VSH  model  is  also  compared  with  output  from  the  Mass  Spectrometer- 
Incoherent  Scatter  (MSIS-S6)  (Hedin,  1987)  empirical  model. 

The  effects  of  wind  drag  on  satellites  is  quantified  by  comparing  the 
satellite  drag  due  to  in-track  winds  to  the  drag  without  winds  for  different 
cases  dependent  on  magnetic  activity  and  satellite  orbit  paths. 
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The  first  chapter  contains  background  information  on  the  theory  of 
upper  atmospheric  dynamics  and  the  parameters  that  control 
thermospheric  density  structure.  Also,  increases  in  drag  due  to  polar 
winds  are  explained. 

The  next  chapter  contains  details  on  the  research  tools  utilized. 
Information  on  thermospheric  general  circulation  models  (TGCMs)  is 
included  as  a  prelude  to  the  section  on  the  VSH  computer  model.  VSH 
model  output  is  compared  against  total  density  measurements  from  the 
SETA-2  instrument  and  the  instrument  on  the  Dynamics  Explorer  2  (DE-2) 
satellite.  The  1986  version  of  the  MSIS  empirical  model  (Hedin,  1987)  is 
used  to  make  model-to-model  comparisons.  And  finally,  the  technique 
used  to  compare  satellite  data  to  the  models  is  illustrated. 

The  next  chapter  goes  into  the  model  validation  program  that  was 
created  and  the  systematic  process  used  to  compare  model  and  data. 
Satellite  data  was  binned  using  a  variety  of  parameters  to  evaluate  model 
characteristics.  Inclusion  of  an  additional  solar  flux  parameterization  was 
added  based  on  results  from  this  study.  Validation  studies  such  as  this  are 
important  in  continuously  updating  and  improving  model  performance. 

The  atmospheric  drag  effects  of  wind  on  a  satellite  is  detailed  in 
chapter  four.  And  finally,  a  discussion  of  current  and  future  research  is 
summarized  in  chapter  five. 

1.2  Thermospheric  neutral  density  and  winds 

Since  the  objective  of  this  thesis  is  the  validation  of  neutral  density 
output  of  the  VSH  thermospheric  model,  some  background  on  the  modelled 
region  should  be  provided. 
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In  general,  the  mean  atmospheric  density  decreases  with  altitude. 
In  the  lower  atmosphere  up  to  100  km  most  of  Jie  chemical  species  behave 
as  one  and  decrease  exponentially  by  the  barostatic  distribution  equation 

p(z)  =  p(zQ)e  *s)  (1.1) 

where  density,  p,  drops  off  with  respect  to  some  reference  altitude  z0 
exponentially  by  the  scale  height  of  the  gas.  The  scale  height 

IcT 

H  =  ~  (1.2) 

mg 

of  the  overall  gas  with  mean  molecular  mass  in  at  a  constant  temperature, 
T,  and  constant  gravity,  g,  (k  being  Boltzmann’s  constant)  is  the  “thickness” 
that  the  gas  would  exhibit  if  held  at  a  reference  pressure.  Based  on  this 
equation  it  can  bo  seen  that  h  savy  gases  have  small  scale  heights  and  light 
gases  b*ve  large  scale  heigh .‘s.  This  lower  region  is  called  the  homosphere 
and  is  ^emanated  by  eddy  diffusion,  the  process  by  which  the  individual 
gases  are  fully  mixed  due  to  winds  and  turbulence. 

Above  100  km,  well  into  the  thermosphere,  eddy  diffusion  gives  way  to 
molecular  diffusion  in  the  heterosphere.  The  individual  species 
gravitationally  separate  based  on  their  own  masses  and  exhibit  individual 
scale  heights.  In  equation  1 .2,  H  and  m  change  to  H,  and  mj  whore  the 

subscript  i  represents  each  individual  species  scale  height  and  mass, 
respectively.  Figure  1.1  illustrates  this  distinction  in  scale  heights  in  the 
homosphere  and  heterosphere.  It  can  be  seen  that  the  heavier  species  such 
as  argon  (Ar)  and  ozone  (O3)  begin  to  fall  off  (with  respect  to  density)  more 

rapidly  than  the  lighter  species  (He)  around  100  km.  In  the  heterosphere 


figtiro  i.i  Density  pfofttec  ol  various  atmospheric  constituents. 


most  of  the  species  are  in  molecular  equilibrium,  meaning  that  there  is  no 
net  motion  of  an  individual  species  through  tho  rest  of  the  gases. 

The  most  important  and  abundant  chemical  species  in  the 
thermosphere  are  O,  O2,  and  N2.  Molecular  nitrogen  is  chemically  inert 

and  is  in  molecular  equilibrium,  and  therefore,  is  controlled  by  dynamic 
transport  However,  0  and  O2  are  not  in  molecular  equilibrium,  instead, 

their  density  profiles  are  controlled  by  chemistry  in  the  lower 
thermosphere. 


L2J  Control  of  thermospheric  density 
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The  density  structure  of  the  thermosphere  is  best  described  by 
evaluating  one  of  the  primitive,  or  governing,  equations  for  the 
thermosphere.  There  are  three  equations,  the  conservation  of  mass, 
momentum,  and  energy.  The  equation  for  the  conservation  of  mass,  also 
called  the  continuity  equation,  best  describes  the  processes  responsible  for 
the  composition  and  density  structure  of  the  thermosphere. 

The  single  species  continuity  equation  for  the  thermosphere  is  given 
in  terms  of  the  ratio  of  the  individual  species  density  to  that  of  the  whole 
gas,  the  mass  mixing  ratio  SP: 


d'F 

at 


£1A 

t  3z 


m  fTpp 


+  e 


dz 


K(z)e 


d  v7 
dz 


T  -  V 


dz 


+S-R 


(1.3) 


The  terms  on  the  right-hand  side  of  the  equation  are,  respectively,  changes 
in  comnocition  due  to  molecular  diffusion,  eddy  diffusion,  horizontal  and 
vertical  advection,  and  finaUy,  chemical  production  and  loss  (Killeen  et  al., 
Iy88).  The  various  terms  are: 

V  =  horizontal  velocity  vector 
w  a  vertical  neutral  velocity 
T  -=  temperature 

T00=  s  u  if  ace  temperature  (273  K) 
x  =  diffusion  time  scale 
L  =  matrix  of  diffusive  f'^uilibrium  solutions 
a  «  matrix  of  diffusion  coe^icionts 
S  =  matrix  for  chemical  production  rates 
R  -  matrix  for  chemical  loss  rate^ 
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Looking  at  each  term  and  applying  it  to  the  lower  thermosphere,  we 
can  choose  the  most  important  process  that  affects  variations  in  the  total 
density. 

The  chemical  sources  and  sinks  are  important  for  individual  species 
chemistry  but  not  so  for  total  density.  Also,  molecular  and  eddy  diffusion 
have  significance  for  individual  species  density  profiles  but  are  of  little 
importance  in  specifying  perturbations  in  total  density  (Killeen  et  al.,  1988). 
Therefore,  variations  in  the  mean  density  structure  of  the  thermosphere 
are  primarily  controlled  by  horizontal  and  vertical  advective  effects  (Killeen 
et  al.,  1988).  These  advective  effects  are  caused  by  heating  or  momentum 
exchange. 

Thermospheric  neutral  winds  respond  to  differential  heating  by 
flowing  from  the  warm  day  side  to  the  cooler  night  side  of  the  earth  (figure 
1.2)  and  from  the  summer  to  the  winter  hemisphere.  This  thermospheric 
flow  results  in  restoring  the  global  temperature  distribution.  Subsequently, 
this  flow  causes  day-night  and  summer-winter  horizontal  density 
gradients  due  to  the  advective  effects  of  the  winds. 

Winds  also  reside  in  large  scale  systems  such  as  planetary  waves 
and  diurnal  tides  (the  “breathing  effect"  of  the  atmosphere  in  response  to 
the  daily  tolar  EUV/UV  heating),  and  small  scale  systems  such  as  gravity 
waves  and  localized  upwelling  (as  a  result  of  local  heating)  (Killeen  et  at., 
1988).  These  occurrences  of  momentum  transfer  create  vertical  winds, 
resulting  in  changes  to  the  density  profile.  For  instance,  the  diurnal  tide 
creates  upward  winds  (increased  total  mass  density)  during  the  day  and 
downward  winds  (decreased  total  mass  density)  during  the  night.  Figure 
1.3  clearly  illustrates  this  diurnal  variation,  where  the  vertical  axis  is 


,  £■ 


figure  1.2  Illustration  of  horizontal  wind  flow  due  to  differential  heating 
(Hedin  et  al..  1988). 
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departure  from  the  U.S.  standard  atmosphere  density  and  the  horizontal 
axis  is  local  time.  The  thermospheric  wind  system,  and  resultingly,  the 
density  structure,  can  also  be  disrupted  by  the  11  year  solar  cycle  and 
geomagnetic  storms.  More  information  on  geomagnetic  storms  and  their 
effect  on  total  mass  density  is  discussed  in  section  1.2.1 .2.  These  are 
specific  occurrences  which  essentially  result  in  changes  to  the  density 
structure  described  in  the  above  paragraphs.  These  features  are 
perturbations  that  cause  variations  in  the  mean  density  structure,  but  what 
restores  the  mean  density  structure? 

For  a  “steady  state”  thermosphere,  the  horizontal  and  vertical  total 
mass  density  profiles  are  governed  by  the  distribution  law,  equation  1 .1 ,  for 
each  of  its  individual  species,  which  simply  says  that  the  temperature 
profile  is  the  controlling  factor.  The  temperature  profile,  in  turn,  is 
controlled  by  thermospheric  heat  sources  and  sinks  which  are  solar 
EUV/UV  heating  (source)  and  downward  heat  conduction  (sink). 
Downward  heat  conduction  to  the  mesosphere  is  the  primary  means  of  the 
thermosphere  to  dissipate  heat.  Therefore,  the  balance  between  these  two 
factors  control  the  mean,  or  globally  averaged,  density  structure. 

L2X1  Solar  flux  and  geomagnetic  indices 

The  10.7  cm  (2800  MHz)  solar  radio  emission  has  been  extensively 
used  to  indicate  solar  activity.  This  emission  has  been  designated  the  FI 0.7 
index  and  correlates  well  with  the  sunspot  number  (figure  1.4),  which, 
based  on  hundreds  of  years  of  recorded  observations,  is  clearly  indicative  of 
solar  activity  (AFGL,  1986). 


MEAN  SU 


An  FI  0.7  value  gives  the  day-to-day  variability,  and  can  vary  greatly 
over  a  short  period.  A  slowly  variable  index  was  found  to  be  useful  in  upper 
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figure  1.4  Illustration  showing  correlation  between  sunspot  numbers  to  solar  flu* 
(AFGL,  1985). 
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atmospheric  models;  the  F10.7  average  (F10.7a)  value,  accounts  for  the 
long-term  variability  (Hinteregger,  1981)  of  the  sun.  This  value  is  a 
running  81  day  average  of  the  previous  FI 0.7  values  (Hinteregger,  1981). 
For  some  thermospheric  computer  models,  as  in  the  case  of  the  MSIS 
model,  the  FI 0.7a  value  is  calculated  by  a  ±40  day  mean  value  for  its 
determination  of  exospheric  temperatures  (Hinteregger,  1981).  A  study  by 
Hinteregger  [1981]  using  AE-E  observations  lead  to  the  conclusion  that  EUV 
fluxes  correlate  with  FI 0.7a  more  so  than  with  the  daily  F10.7  value.  This 
discovery  also  relates  well  with  the  results  of  this  validation  study,  showing 
that  VSH  needed  Fl0.7a  as  an  input  parameter  in  addition  to  the  daily  F10.7 
value  to  improve  its  density  output. 

Even  though  the  10.7  cm  solar  flux  is  a  good  indicator  of  the  sun’s 
variability,  it  still  can  only  be  considered  a  proxy  to  the  more  indicative  solar 
EUVAJV  emissions.  Tho  10.7  cm  emission  is  different  from  these  stronger 
UV  emissions.  While  the  former  comes  from  the  upper  chromosphere, 
lower  corona  region  of  the  sun,  EUVAJV  emissions  come  from  a  broad 
region  of  the  solar  atmosphere  where  the  temperatures,  densities,  and 
altitudes  are  quite  different  (AFGL,  1985). 

EUVAJV  emissions  are  primarily  responsible  for  the  earth’s  global 
energy  budget  (Torr  et  al.,  1981)  and  would  be  the  perfect  indicator  for  solar- 
terrestrial  interactions.  However,  trying  to  index  these  many  different 
EUVAJV  spectral  lines  would  make  thoir  use  impractical.  Also,  the 
success  of  the  MSIS  model  and  Jacchia  satellite  drag  model  and  their  lovel 
of  verified  inconsistencies  with  actual  aeronomical  observations  shows  that 
a  large  number  of  solar  variability  indices  is  not  really  needed  anyway 
(Hinteregger,  1981). 
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A  number  of  indices  have  been  developed  to  measure  geomagnetic 
activity.  The  K  index  was  developed  to  show  irregular  variations  in 
magneto  grams,  which  are  indicators  of  magnetic  disturbances  such  as 
geomagnetic  storms. 

The  K  index  is  taken  at  3  hour  intervals  based  on  the  difference 
between  the  highest  and  lowest  deviations  from  the  regular  daily  variation 
in  the  magnetic  field.  Due  to  the  structure  of  the  earth’s  magnetosphere, 
this  index  primarily  indicates  auroral,  'or  high-latitude,  activity  (AFGL, 
1985). 

K  values  are  roughly  logarithmic  integers  from  0  to  9.  The  Kp  index, 
where  p  stands  for  “planetary",  is  intended  to  reflect  the  worldwide  average 
level  of  magnetic  activity  (AFGL,  1985). 

The  ap  and  Ap  index  is  also  widely  used,  mostly  when  averaging  of 
the  3-hourly  Kp  values  is  needed  to  get  a  daily  average.  The  ap  value  has  a 
linear  range  that  is  converted  from  the  Kp  range  by  use  of  a  standardized 
table  (AFGL,  1986).  The  ap  values  are  then  averaged  over  eight  3-hour 
intervals  to  get  the  daily  Ap  index.  Ap  values  range  from  0  to  400. 

The  need  for  the  indexing  of  solar  flux  and  geomagnetic  activity 
primarily  arises  in  computer  models  that  simulate  the  earth’s  upper 
atmosphere.  The  VSH  model,  for  instance,  is  strongly  dependent  on 
variations  in  the  daily  solar  flux,  and  a  simple  index  characterizing  this 
field  helps  in  simplifying  the  model. 

L2X2  Geomagnetic  storms 

Geomagnetic  disturbances  come  about  from  the  transfer  of  energy 
from  the  solar  wind  to  the  earth’s  magnetosphere.  The  largest 
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disturbances  of  the  magnetosphere  are  called  magnetospheric  storms  and 
the  resulting  disturbance  of  the  geomagnetic  field  is  called  a  geomagnetic 
storm  (AFGL,  1985).  These  larger  disturbances  are  primarily  due  to  a 
sudden  impulse  from  a  shock  in  the  solar  wind.  Shocks  in  the  solar  wind 
are  sometimes  caused  by  solar  flares  (Roble  et  al,  1987). 

This  transfer  of  energy  causes  enhanced  particle  precipitation  (ie. 
enhanced  mean  particle  energy)  in  the  earth’s  auroral,  or  high  latitude, 
region  and  also  causes  magnetospheric  convection.  These  processes 
disturb  the  high-latitude  thermospheric  wind  circulation  patterns, 
temperature,  and  composition. 

Forbes  et  al.  [1987]  describes  the  effects  on  total  mass  density  by  such 
a  storm  using  satellite  neutral  wind  and  total  mass  density  measurements. 
Their  results  show  that  densities  in  the  high  latitude  region  of  the  northern 
hemisphere  were  enhanced  on  the  day-side  by  about  60%  during  the  peak  of 
the  magnetic  storm.  This  enhancement  corresponded  with  intense 
heating.  The  measurement  was  taken  at  180  km  at  the  throat  of  the  ion 
convection  pattern  (see  section  1.2.2).  This  area  of  the  ion  convection 
pattern  has  enhanced  particle  and  Joulo  heating  and  strong  convergence 
during  geomagnetic  storms  (Roble  et  al.,  1987).  Joule  heating  is  the 
frictional  heating  due  to  the  changes  between  the  accelerated  ions  and  the 
slower  neutral  particles  that  flow  with  them. 

The  disturbance  increases  the  density  as  it  propagates  from  tine  high 
latitudes  to  the  equator.  Satellite  measurements  on  the  night-side  show 
less  latitudinal  density  enhancement,  where  the  disturbance  did  not 
propagate  equatorward  as  extensively  (Forbes  et  al.,  1987).  After  the  storm 
passes,  densities  continue  to  stay  high  (10-20%)  (Forbes  et  al.,  1987) 
compared  to  pre-storm  conditions. 
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XSL2  High  latitude  winds 

The  previous  section  describes  thermospheric  composition  and  winds 
in  general.  I  will  now  go  into  some  background  information  needed  to  help 
explain  my  research  in  the  effects  of  atmospheric  drag  due  to  winds 
experienced  by  satellites. 

Thermospheric  winds  reach  a  maximum  in  the  polar  cap  regions  of 
the  earth.  During  strong  geomagnetic  storms,  neutral  winds  can  roach 
speeds  well  over  1000  m/sec  (Killeen  et  al.,  1983),  and  when  compared  with 
the  average  satellite  speed  of  8000  m/sec,  the  wind  can  change  the  drag/lift 
forces  on  the  satellite  by  50%  under  worse  case  conditions. 

In  the  higher  latitudes  neutral  winds  follow  the  ion  convection 
pattern.  This  circulation  pattern  results  from  interactions  between  the 
solar  wind,  magnetosphere  and  the  ionosphere. 

The  solar  wind  interacts  with  the  magnetosphere  producing  dynamo 
electric  fields.  An  external  dynamo  is  produced  that  creates  dawn  to  dusk 
electric  fields  (EpC)  that  map  down  to  the  polar  cap  region.  An  internal 
dynamo  is  also  produced  that  creates  electric  fields  (Ea)  that,  when  mapped 
down  to  the  polar  cap,  are  reversed  in  direction,  creating  a  dusk  to  dawn 
electric  field.  Figure  1.5  (Thayer,  1990)  shows  the  location  and  direction  of 
EpC  and  Ea  and  the  resulting  wind  pattern  (black  arrows)  that  develops 

around  these  electric  fields. 

These  electric  fields  force  ions  into  motion  creating  this  two-cell  ion 
convection  pattern  that  flows  anti-sunward  in  the  polar  cap  with  a  sunward 
return  flow  on  the  outer  boundaries.  Due  to  ion-neutral  collisional  drag,  or 
simply  called  ion  drag,  the  neutrals  are  "dragged’’  along  following  this 
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figure  t.S  Schematic  of  the  resulting  circulation  cells  from  externally  imposed 
electric  fields  (Thayer,  1990}. 

pattern.  Neutral  wind  wind  speeds  in  this  pattern  are  normally  100-300 
m/sec  across  the  polar  cap. 

This  pattern  is  generally  oriented  along  the  0000-1200  local  time 
meridian  (Marcos,  1985).  But  during  strong  geomagnetic  storms  the 
pattern  expands  to  lower  latitudes,  becomes  more  antisymmetric  and  can 
change  its  orientation.  Neutral  winds  can  also  increase  during  these 
storms  to  over  1000  m/s ec.  The  enhar  "nd  circulation  pattern  prevails  for 
hours  after  a  storm  ceases  (Marcos,  1985). 


1.3  Aerodynamic  drag 
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The  total  mass  density  and  its  variability  is  directly  related  to  the 
atmospheric  drag  effect  on  satellites.  Therefore,  our  ability  to  measure  or 
model  the  mean  density  and  its  variability  is  important  to  the  success  of  any 
satellite  mission.  Atmospheric  drag  should  be  taken  into  consideration  for 
orbital  lifetime,  design  and  control,  tracking  and  on-board  fuel 
requirement  (Marcos  et  al.,1989).  Such  empirical  models  as  Jacehia-70  and 
MSIS  have  been  and  still  are  used  successfully  to  estimate  the  impact  of 
drag  on  satellites.  However,  accuracies  of  these  models  have  remained  at 
~15%  for  the  past  two  decades  (Marcos,  1991).  Soon,  the  VSH  model  will 
also  be  used  for  this  purpose,  but  with  improved  accuracy  and  possibly  in  a 
forecast  mode  of  operation. 

Although  there  has  been  much  research  into  the  aerodynamic  drag 
effect  on  satellites  due  to  atmospheric  density,  there  has  been  little  study  of 
the  effect  that  winds  by  themselves  have  on  satellite  drag.  Noutral  winds 
directed  opposite  to  that  of  the  satellite  velocity  vector  have  the  strongest 
drag  effect  on  satellites,  since  the  drag  force  is  proportional  to  the  square  of 
the  relative  velocity. 

Studies  by  Marcos  (1985)  have  shown  that  the  drag  force  on  a  satellite 
can  increase  by  -5%  for  every  200m/soc.  This  applies  only  if  the  winds  are 
anti-parallel  to  the  direction  of  the  satellite.  During  geomagnetic  storms, 
enhanced  winds  in  the  twin  circulation  cells  can  increase  the  drag  on  a 
satellite  significantly. 

I  attempted  to  expand  upon  this  study  by  determining  the  integrated 
drag/lifl  forces  (along  the  orbit  track)  on  a  satellite  due  to  neutral  winds.  I 
used  a  typical  DE-2  orbit  and  computed  neutral  winds  along  the  satellite 
track  with  the  VSH  model.  It  was  found  that  742  m/sec,  anti-parallel  winds 
changed  the  satellite  drag  by  23%. 
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The  equation  used  to  determine  aerodynamic  drag  is  similar  to  that 
used  in  section  2.4.2  but  with  a  few  assumptions: 

Fn  =  pV2  (1.4) 

The  vector  FD  represents  u.e  drag  force  and  p,  the  total  mass  density.  V  is 
the  atmospheric  mass  velocity, Vw,  relative  to  the  satellite  velocity,  V,: 

V  =  V.-VW  (1.4) 

This  is  different  from  the  expression  in  section  2.4.2  because  we  are 
making  relative  comparisons,  so  to  simplify  the  parameterization  in  the 
FORTBAN  program  we  assume  all  the  other  parameters  (satellite  frontal 
area  and  drag  coefficient)  are  simply  1.  The  factor  of  1/2  is  also  taken  out 
for  simplification. 

I  attempted  to  do  more  research  in  this  matter,  calculating  the 
integrated  drag  along  the  satellite  orbit,  and  determining  exactly  the 
contributions  to  drag  in  the  polar  cap  region.  However,  due  to  time 
constraints  I  was  required  to  cut  the  study  short,  providing  only 
information  on  the  increase  in  drag  at  particular  points  along  the  orbit. 
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CHAPTER 2 


Di 


S3PI1QN  OF  TOOLS  USED 


The  following  chapter  provides  background  information  on  the 
dovices  used  to  conduct  the  VSH  validation  study.  A  short  history  on 
thermospheric  general  circulation  models  (TGCMs)  is  included  since  the 
VSH  model  is  based  on  the  output  of  the  NCAR-TIGCM.  And,  of  course, 
details  on  the  techniques  used  in  the  VSH  model  are  important  in 
understanding  the  nature  of  the  computer  model’s  output.  The  MSIS 
empirical  model  is  referenced  by  both  TGCM  and  VSH  models  under 
various  situations  and  is  also  used  in  this  study  comparing  results  to  the 
SETA  measurements.  The  SETA  instrument  and  it’s  measurements  used 
in  the  validation  procedure  are  explained  also.  DE-2  density  measurements 
were  also  used  to  a  small  extent  to  possibly  provide  altitudinal  comparisons 
of  the  models  in  the  study.  So  background  r>n  the  neutral  atmosphere 
composition  spectrometer  instrument  is  also  included. 


2,1  Thermospheric  General  Circulation  Models  (TGCMs) 


Historically,  development  of  general  circulation  models  has  been 
limited  due  by  the  inability  of  computers  to  solve  the  primitive  equations 
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associated  with  aeronomieal  applications  in  a  timely,  usable  manner. 
After  years  of  developing  and  using  one  dimensional,  then  two  dimensional 
global  thermospheric  models,  it  was  finally  realized  that  to  best  simulate 
the  thermosphere,  a  self-consistent  three  dimensional  model  that 
simultaneously  solved  the  applicable  conservation  equations  was  needed. 

The  National  Center  for  Atmospheric  Research  Thermospheric 
General  Circulation  Model  (NCAR-TGCM)  (Dickinson  et  al.  1981)  and  the 
University  College  at  London  Thermospheric  General  Circulation  Model 
(Fuller-Roweil  and  Rees,  1980)  were  the  first  models  to  solve  the  three 
dimensional  time  dependent  equations  that  used  to  describe  the 
thermosphere.  The  NCAR-TGCM  is  based  on  a  circulation  model  of  the 
lower  atmosphere,  used  in  forecasting  tropospheric  weather,  with  added 
physical  processes  applicable  to  the  upper  atmosphere  (Dickinson  et  al., 
1981).  The  development  and  use  of  global  circulation  models  like  these  has 
helped  provide  insight  into  the  dynamics  of  the  upper  atmosphere. 

The  resolution  of  the  NCAR-TGCM  is  5°  in  latitude  and  longitude 
and  has  24  pressure  levels  which  covers  altitudes  from  about  90  km  to  500 
km.  The  TGCM  employs  empirical  and  semi-empirical  models  and 
parameterizations  to  provide  its  output  (Dickinson  et  al,  1981).  Output 
includes  neutral  winds,  temperature,  and  composition. 

The  dominant  forcing  of  thermospheric  processes  is  solar  radiation 
in  tlie  extreme  ultra-violet  (EUV)  and  ultra-violet  (UV)  range.  Dickinson  et 
al.  [1981]  describes  using  solar  EUV  fluxes  from  Hinteregger  [1981]  and  UV 
fluxes  form  Torr  et  al.  [1980]  to  characterize  the  solar  heating  and 
photodissociation  terms.  To  represent  the  electron  density  distribution 
from  90  to  500  km,  the  ionospheric  model  by  Chiu  [1975]  was  used 
(Dickinson  et  al.,1981).  Neutral  composition  and  temperature  to  calculate 
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the  global  average  is  taken  from  the  MSIS  empirical  model  from  Hedin 
[1977]  (Dickinson  et  al.,  1981). 

In  recent  years  the  original  NCAR-TGCM  has  become  more  self- 
consistent,  relying  less  on  empirical  models,  and  broadened  it’s 
boundaries.  A  Thermosphere/Ionosphere  General  Circulation  Model 
(TIGCM)  has  been  developed  that  couples  the  dynamics  between  the  two 
regions,  alleviating  part  of  the  need  for  the  Chiu  [1975]  model.  Included  in 
a  recent  version  of  the  model  are  the  effects  of  neutral  constituents  on  ions 
and  has  led  to  the  development  of  the  Thermosphere/Ionosphere 
Electrodynamics  General  Circulation  Model  (TIEGCM).  Finally,  the 
extension  of  the  model  to  include  the  chemistry  and  dynamics  of  the 
mesosphere  will  lead  to  the  TIMEGCM  (communication  with  R.  Roble, 
1991). 

2L2  Vector  Spherical  Harmonic  (VSH)  Model 

The  Vector  Spherical  Harmonic  (VSH)  semi-empirical  model  is 
based  on  output  from  the  NCAR-TIGCM.  Input  needed  to  drive  the  VSH 
model  include  temporal,  spatial  and  geophysical  parameters  (Killeen  et  al., 
1991),  The  model  was  developed  partly  to  give  researchers  a  simple, 
relatively  easy  to  use  computer  model  of  the  global  thermospheric  structure 
that  can  run  on  a  small  computer.  The  basic  idea  of  the  VSH  model  is  to 
take  a  spectral  expansion  in  latitude  and  longitude  of  the  gridded  output 
from  the  NCAR-TIGCM  using  scalar  and  vector  spherical  harmonics 
(Thayer,  1990).  Its  use  of  vector  spherical  harmonics  if  where  the  model 
gets  its  name. 
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The  specification  of  neutral  density  is  very  important  in  determining 
the  atmospheric  drag  on  artificial  satellites.  Historically,  the  accuracy  of 
the  prediction  of  neutral  density  has  been  limited  to  no  better  than  15%.  The 
VSH  model's  goal  is  to  reduce  this  error  to  5%  by  using  real-time  satellite 
data  to  refine  the  model’s  output  (Killeen  et  al.,  1991). 

VSH  is  referred  to  as  a  semi-empirical  model  because  it  employs 
output  from  the  MSIS-86  empirical  model  for  altitudes  above  500  km  and 
altitudes  below  90  km. 

The  TIGCM  output  fields  (winds,  temperature  and  composition)  are 
in  5  degree  grids  in  latitude  and  longitude  over  24  pressure  levels  vertically 
at  selected  Universal  Times  (UTs).  A  procedure  is  conducted  on  the 
TIGCM  history  file  which  performs  the  spectral  expansion  described  above 
using  vector  spherical  harmonics  on  the  wind  fields  and  scalar  spherical 
harmonics  on  the  temperature  and  composition  output  fields  of  the 
TIGCM.  The  coefficients,  or  "weighting  factors",  created  due  to  this 
expansion  are  placed  in  a  coefficient  file.  The  coefficients  are  used  by  the 
VSH  model,  by  various  interpolation  methods,  to  fit  smooth  curves  over  the 
TIGCM  output  fields  to  show  variation  in  winds,  temperature  and 
composition  with  respect  to  the  initial  conditions.  Inputs  include  temporal 
(UT),  spatial  (latitude,  longitude)  and  geophysical  conditions  (Ap,  solar 
flux)  (Killeen  et  al.,  1991).  This  procedure  results  in  converting  discrete 
gridded  TIGCM  output  to  continuous  output  used  by  VSH  (Killeen  et  al., 
1991). 

The  coefficient  files  mentioned  above  are  placed  in  a  coefficient 
library.  The  coefficient  library  is  produced  from  a  number  of  runs  of  the 
NCAK-TIGCM  and  are  limited  to  keep  the  storage  requirements  at  a 
minimum.  These  runs  are  based  on  a  broad  range  of  geophysical 
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conditions  including  seasonal  variation,  magnetic  activity,  and  solar 
activity.  A  matrix  of  runs  exist  that  cover  three  time  periods  over  the  year 
(June  and  December  solstice  and  equinox  conditions),  three  levels  of 
magnetic  activity  (Ap  values  of  5, 11  and  32),  and  two  levels  of  solar  activity 
(low  and  high  values  of  FI 0.7a  and  FI  0.7  to  represent  solar  minimum  and 
maximum)  (Killeen  et  al.,  1991),  so  there  are  presently  18  TIGCM  runs  in 
the  coefficient  library  for  use  by  VSH. 

2.2.1  Spherical  Harmonics 

The  following  paragraphs  describe  the  scalar  and  vector  spherical 
harmonics  used  by  the  VSH  model  to  convert  the  TIGCM  history  file  to  the 
coefficient  files  used  by  VSH  to  resynthesize  output  fields.  More  details  on 
the  mathematics  of  spherical  harmonics  can  be  found  in  Thayer  [1990]. 

Since  the  earth  is  primarily  spherical,  one  of  the  best  ways  to  describe 
horizontal  variations  over  the  earth  is  by  using  spherical  harmonics.  This 
method  helps  alleviate  the  problem  of  discontinuities  occurring  over  the 
poles  as  lines  of  longitude  converge.  A  scalar  quantity  (e.g.  density  or 
temperature)  can  be  described  in  scalar  spherical  harmonics  at  a 
particular  location  on  a  sphere  from  the  sum 

q(M)  =  (2.!) 

m.n 

where  6  is  colatitude,  X  is  longitude,  Ynm  is  the  scalar  spherical  harmonic 
and  am  n  are  the  complex  coefficients.  The  value  Y"  can  be  expressed  as: 


Y™(0,A)  =  P„w(cos0)  eimX 


(2.2) 
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with  P®  being  the  associated  Legendre  functions 

1  J(m+n) 

P“m(0)  =  2^!(8ine)m^(xa"1)n;  X  “  0086  (23) 

The  top  index  m  is  the  order  or  zonal  wavenumber  at  a  particular  latitude 
and  the  bottom  index  n  is  the  degree.  The  latitude  wavenumber  is 
determined  by  n-m  (Thayer, 1990). 

The  vector  quantity  V,  in  this  case  the  TIGCM  output  wind  field,  is 
expanded  using  vector  spherical  harmonics  which  has  this  form: 

(2.4) 

m,n 


where  a,b,c  are  the  complex  coefficients,  and 


(2.6  a„b,c) 


where  A®  and  B®  are  related  to  the  Legendre  functions  P0m  by 


m  pm 

sin  9  0 


(2.6  a, b) 
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(Killeen,  1987). 

The  complex  coefficients  that  go  into  the  library  are  computed  by 
taking  a  least  squares  fit  of  the  TIGCM  wind  output  (Killeen  et  al.,  1991). 
Then  by  putting  these  coefficients  back  into  equation  1.4,  the  global  wind 
field  can  be  resynthesized  with  continuous  output  (Thayer,  1990). 

2^2  Interpolation  Methods 

The  curves  described  above  are  not  only  fitted  through  a  single 
TIGCM  run  but  can  also  be  fitted  across  different  TIGCM  runs,  such  as  the 
manner  solar  flux  values  are  interpolated,  to  better  represent  the  initial 
conditions.  The  interpolation  methods  used  by  VSH  are  described  in  the 
following  paragraphs. 

The  interpolation  method  used  to  represent  vertical  structure  is  a 
cubic  spline,  which  is  a  piecewise  polynomial  (Killeen  et  al.,  1991).  The 
vertical  structure  of  the  atmosphere  varies  greatly  due  to  physical  processes 
that  occur  at  different  altitudes.  The  vertical  profiles  of  density  and 
temperature  vary  significantly  between  100  km  and  500  km. 

Near  100  km,  where  the  turbopause  exists,  density  falls  off  in  a 
uniformly  exponential  manner  below  the  turbopause.  Above  the 
turbopause,  the  individual  chemical  species  begin  to  separate  and  their 
density  drops  off  depending  on  their  individual  mass.  Near  this  same 
altitude  at  the  mesopause  the  temperature  of  the  atmosphere  is  at  its 
coldest  and  quickly  increases  to  its  maximum  near  500  km. 

In  the  VSH  model  cubic  splines  are  used  in  the  altitude  range  from 
110  -  500  km.  Above  that  altitude  temperatures  from  the  MSIS  model  are 
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used.  Coefficients  of  the  polynomial  are  different  for  each  subinterval  of  the 
altitude  range  (Killeen  et  al.,  1991  \ 

A  Fourier  time  series  is  used  to  reproduce  diurnal,  semidiurnal,  and 
terdiumal  variations.  A  Fourier  representation  assumes  a  periodic  nature 
of  that  field.  This  is  generally  true,  for  example,  for  diurnal  tides,  but 
during  magnetic  storms,  fields  can  vary  significantly  before  returning  to  a 
normal  state.  Plans  exist  to  implement  the  use  of  a  spectral  fit  for 
magnetic  storm  cases  (Killeen  et  al.,  1991). 

The  Fourier  series  used  is  expressed  as: 


am,n.«  =  ^“+Zak.»^C°8kt+  Zk.m,»..Binkt 
*  k»l 


T 

I 

k.l 


(2.7) 


For  T=3  this  represents  a  constant  term,  three  symmetric  coefficients,  and 
three  antisymmetric  coefficients  (diurnal,  semidiurnal,  and  terdiurnal) 
(Killeen  efcal.,1991). 

To  account  for  seasonal  variations  over  the  year  a  trigonometric 
interpolation  method  is  used.  Solstice  and  equinox  conditions  due  to  the 
annual  solar  cycle  are  represented  by  a  sinusoidal  curve  of  the  radiation 
absorbed  by  the  earth  (Killeen,  et  al.,  1991). 

Previously,  the  VSH  model  used  only  the  FI 0.7  solar  EUV  proxy  in 
representing  solar  flux  variations  (Killeen,!  991).  As  mentioned  before, 
solar  EUV  and  UV  radiation  affects  thermospheric  dynamics  greatly,  so  a 
good  representation  of  solar  flux  is  very  important  to  the  VSH  model. 
During  the  model  validation  process  described  in  chapter  3,  VSH  model 
density  was  noticed  being  greatly  overestimated  under  certain  conditions. 
It  was  discovered  that  this  overestimation  was  occurring  during  times  of 
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strong  solar  flux,  or  high  FI  0.7  values.  A  better  representation  of  solar  flux 
was  needed. 

The  FI  0.7a  value  (the  81  day  running  average)  is  very  important  in 
describing  low  altitude  dynamics,  while  the  daily  solar  flux  value,  F10.7,  is 
more  important  in  the  higher  altitude  region.  A  parameterization 
combining  the  effects  of  these  two  values  into  one  value,  FI 0.7’,  is  used  to 
take  into  account  this  height  dependence 

FI 0.7’  =  FlO^a  +  k^h^CFlO.T-FlO^aJ+kjChMFlOJ-FlO^a)8  (2.8) 

Here  ki  and  k2  are  height  dependent  coefficients  derived  from  the  MSIS 
model.  The  weights  increase  monotonically  with  height  thereby  giving 
more  importance  to  the  FI  0.7  value  at  high  altitudes.  This 
parameterization  corrected  the  density  overestimation  and  resulted  in 
improved  density  output. 

Geomagnetic  activity  is  indexed  by  the  Ap  value.  Since  the  Ap  index 
is  a  linear  measure  of  the  magnetic  activity,  variations  of  magnetic  storm 
activity  are  represented  by  a  linear  interpolation  across  Ap  (Killeen  et  al., 
1991). 


2L2L3  Truncation 

One  of  the  important  precepts  of  the  VSH  model  is  to  retain  the  high 
latitude  resolution  feature  of  the  TIGCM  output  fields  (e.g.  the  vector  wind 
fields)  but  keeping  the  coefficient  library  down  to  a  manageable  size.  This 
is  done  by  truncating  the  coefficients  by  a  certain  amount  which  is 
dependent  on  the  indices  n  and  m.  To  get  good  resolution  from  the 
spherical  harmonics  the  latitudinal  wavenumber  (n-m)  must  be 
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considerably  larger  than  the  zonal  wavenumber  (m)  as  in  describing  the 
latitudinal  variations  of  the  wind  fields  at  high  latitudes  (Thayer,  1990). 
This  makes  sense  physically  since  most  of  the  variability  in  a  wind  field,  for 
instance,  occurs  latitudinally. 


LOGto  AMPLITUDE  VSH  COEFFICIENTS 


figure  2.1  Contours  of  the  logarithmic  amplitude  of  the  vector  spherical  harmonic 
coefficients  with  the  three  different  levels  of  truncation  (Kllieen,  1 987). 


Figure  2.1  describes  the  distribution  of  power  (amplitude)  in  the  VSH 
spectral  coefficients  for  a  steady  state  TGCM  run,  under  the  conditions 
stated  in  the  figure,  where  the  pressure  level  Z=1  corresponds  to  a  height  of 
-350  km  (Killeen  et  al.,  1987).  This  is  a  plot  of  the  order  (m)  versus  the 
degree  (n)  of  the  harmonic.  The  lower  right  section  is  unfilled  since  only 
the  values  n>m  are  nonzero.  It  can  be  seen  that  most  of  the  power  lies  in 
the  lower  wavenumbers  but  is  also  relatively  high  latitudinally.  The  power, 
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or  importance,  of  the  coefficients  decreases  rapidly  with  variation  in 
longitude.  Therefore,  to  have  good  resolution  at  the  high  latitudes,  the  most 
important  coefficients  to  keep  are  with  respect  to  the  degree  n  and  less 
important  with  respect  to  the  order  m.  The  best  level  of  truncation  was 
determined  to  be  for  values  of  n^  25  and  m  <,  5  (Killeen  et  al.,1987),  which  is 
level  B  in  the  figure.  This  maintains  good  resolution,  showing  the  high 
latitude  neutral  wind  patterns,  while  cutting  down  the  size  of  the  coefficient 
files  (Thayer,  1990).  Semi-empirical  models  such  as  MSIS  truncate 
coefficients  at  level  A  (Hedin,  1983).  Level  C  would  result  in  excellent 
representation  of  neutral  wind  fields  but  at  a  cost  of  storage  space  (Killeen 
etal.,  1987). 

When  using  a  rectangular  truncation,  even  after  choosing  optimal 
values  of  n  and  m,  undesirable  coefficients  are  included  due  to  the  nature  of 
the  contours.  To  only  attain  the  coefficients  desired,  the  more  logical 
means  of  truncating  would  be  to  follow  the  path  of  the  desired  contour, 
which  also  results  in  optimizing  required  storage  space.  This  will  be  the 
method  of  truncation  for  future  versions  of  VSH  (private  communication 
with  M.Mirmock). 

&3  Mass  Spectrometer-Incoherent  Scatter  (MSIS)  empirical  model 

There  are  a  large  number  of  models  that  simulate  different  aspects  of 
the  thermosphere.  Most  of  them  are  based  on  data  taken  from  in  situ 
measurements  by  satellites,  rocket  borne  instruments,  or  ground  based 
instruments  such  as  incoherent  scatter  radars.  These  empirical  models 
have  been  very  good  at  organizing  large  data  bases  to  be  used  in  theoretical 
work  and  comparing  other  data.  One  of  the  first  of  such  models  was 
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developed  by  Jacchia  [1965]  which  inferred  total  densities  from  satellite 
drag.  Other  models  were  developed  using  temperature  data  from 
incoherent  scatter  radars  (Carru  et  al.,  1967).  But  no  model  had  been 
produced  that  provided  a  combination  of  thermospheric  parameters  until 
MSIS-77  (Hedin,  1977). 

The  MSIS  empirical  model  is  arguably  the  most  widely  used  model 
for  thermospheric  research.  It  uses  data  taken  from  all  the  methods 
mentioned  above.  The  version  used  in  this  study,  MSIS-86,  uses  data  taken 
from  eight  satellites,  numerous  rocket  probes,  and  five  ground-based 
incoherent  scatter  radar  sites  (Hedin,  1987). 

The  MSIS  model  currently  provides  information  on  temperature, 
density,  and  composition  of  the  atmosphere  from  the  upper  mesosphere  to 
the  upper  thermosphere.  This  information  is  also  provided  for  a  large 
range  of  solar  activity  and  magnetic  activity,  for  seasonal,  geographical, 
and  UT  effects. 

Data  is  sorted  into  bins  based  on  latitude,  local  time,  day  of  year,  UT, 
Ap,  F10.7,  and  altitude.  It  is  then  expanded  using  spherical  harmonics  to 
tie  in  all  the  different  data  sots  by  spatial,  temporal,  and  geophysical 
conditions  (Hedin,  1987). 

A  Bates  model  temperature  profile  is  used  for  the  upper 
thermosphere.  An  inverse  polynomial  (Hedin,  1987)  is  used  in  the  lower 
thermosphere  since  this  is  the  region  where  the  atmosphere  changes  from 
mixing  to  diffusive  equilibrium  and  a  simple  profile  does  not  model  Unis 
area  well.  These  methods  are  used  to  determine  the  density  profile  with 
respect  to  geographical  and  solar/  geomagnetic  parameters  (Hedin,  1987). 

The  expansion  mentioned  above  produces  coefTicients  that  arc  used  to 
fit  curves  over  different  data  sets.  In  Hedin  [1977]  the  coefficients  are 
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determined  by  the  120  km  temperature,  120  km  N2  density,  and  the 
exospheric  temperature  for  both  quiet  and  active  magnetic  cases. 

MSIS  has  coverage  over  a  large  range  of  solar  activity,  with  FI 0.7 
values  between  70  and  200  (Hedin,  1983).  Representation  over  magnetic 
activity  has  also  improved  by  using  the  ap,  or  3  hourly,  geomagnetic  index 
as  opposed  to  using  the  daily  Ap  value  (Hedin,  1983). 

With  the  inclusion  of  temperature  and  composition  data  from  the 
Dynamics  Explorer  satellite,  MSIS  has  improved  its  representation  of  the 
high  latitudes  for  solar  maximum  conditions  (Hedin,  1987).  However, 
calibration  problems  in  the  instrument  on  the  DE-2  satellite  that  measured 
total  density  resulted  in  making  the  data  unreliable  and  a  fitting  parameter 
was  required  (Hedin,  1987). 

Such  models  as  MSIS  are  good  at  accurately  assessing  density 
during  quiet  geomagnetic  activity  at  low  latitude  regions  (with  errors  -12- 
15%)  but  during  geomagnetic  storms  at  high  latitudes,  uncertainty  in 
densities  reach  -50%  (Killeen  et  a). ,1988).  The  above  12-15%  accuracies 
currently  exist  as  a  limit  to  the  ability  of  such  models.  Even  additional  data 
sets  such  as  the  use  of  density  and  temperature  measurements  from  the 
Dynamics  Explorer  docs  not  improve  accuracy  (Killeen  et  ah,  1988).  This 
boundary  is  duo  to  the  inherent  accuracy  of  the  sensors  used  to  gather  data, 
possible  errors  in  the  analytical  procedures  used  to  interpolate  data  and 
limited  data  coverage  (Killeen  et  ah,  1988). 

Another,  more  important,  limitation  to  the  MSIS  model  in  particular 
is  the  temporal  and  spatial  resolution  of  this  model.  The  temporal 
resolution  is  limited  by  the  use  of  the  3-hourly  ap  geomagnetic  index 
(Hedin,  1987).  It  is  known  that  magnetic  disturbances  can  occur  well 
within  this  time  frame  (Killeen,  et  ah,  1988).  Also,  the  spatial  resolution  is 
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limited  by  the  level  of  harmonics  used  by  such  empirical  models  (Killeen  et 
alM  1987). 

As  mentioned  in  section  2.2.3,  level  A  in  figure  2.2,  is  the  typical  level 
of  truncation  used  bv  empirical  models  such  as  MSIS  (Thayer,  1990).  This 
level  controls  the  extent  of  the  resolving  capability  of  the  model  and 
subsequently,  the  accuracy  the  model  may  attain.  This  limitation  in 
accuracy  primarily  affects  the  representation  of  the  high  latitude  region 
where  there  is  a  large  amount  of  structure  (ie.  density  enhancements) 
occurring  over  small  distances. 

Until  the  Dynamics  Explorer  2  satellite  provided  a  large  data  base  of 
global  thermospheric  winds,  no  empirical  model  for  such  winds  had  been 
created  (Hedin  et  al.,  1988).  Previous  wind  data  from  ground-based  optical 
instruments,  incoherent  scatter  radar  sites,  and  rocket  chemical  releases 
did  not  provide  adequate  global  coverage.  In  addition  to  DE-2  data,  winds 
from  satellite-based  accelerometers  (Marcos  and  Forbes,  1985)  have  recently 
been  available. 

As  a  result  of  the  availability  of  this  recent  data  an  empirical 
horizontal  wind  model  (HWM-87)  has  been  developed  (Hedin,  1988)  that 
provides  global  winds  for  the  upper  thermosphere.  This  model  uses  wind 
data  from  two  satellites,  and  is  analogous  to  the  MSIS  mode!  (Hedin,  1988). 
The  DE-2  provided  in-situ  zonal  wind  measurements  from  the  wind  and 
temperature  spectrometer  (WATS)  and  remotely  sensed  meridional  winds 
from  the  Fabry-Perot  Interferometer  (FPI)  (Hedin  et  al.,  1988).  The 
Atmosphere  Explorer  E  (AF-E)  satellites  also  supplied  meridional  and 
zonal  winds  from  the  NATE  instrument  (Hedin  et  al.,  1980). 

Due  to  the  data  available,  the  model  is  limited  to  altitudes  above  220 
km,  and  contains  no  F10.7  effects.  However,  effects  due  to  magnetic 
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variability  are  available.  Similar  to  MSIS,  data  was  sorted  based  on 
latitude,  local  time,  day  of  year,  UT,  Ap,  and  FI  0.7  (Hedin  et  al.,  1988).  An 
expansion  using  vector  spherical  harmonics  creates  coefficients  for  the 
spatial  representation.  No  altitude  dependence  is  included  due  to  the 
limited  data  and  theory  suggests  (Rishbeth,  1972)  that  viscosity  in  the  upper 
thermosphere  reduces  vertical  wind  shears,  so  a  vertical  wind  profile  is  not 
as  important. 

The  HWM  only  roughly  depicts  the  neutral  wind  circulation  pattern 
in  the  polar  region,  but  low-  and  mid-latitude  winds  compare  well  with 
ground-based  measurements  (Hedin  et  al.,  1988). 

2,4  SETA-2  instrument 

The  data  set  used  in  this  validation  study  is  composed  of  in-situ  totrl 
mass  density  measurements  from  the  Satellite  Electrostatic  Triaxial 
Accelerometer  (SETA)  instrument.  It  was  flown  on  a  low  altitude  satellite 
with  an  83  degree  inclination  during  1 982.  This  instrument  has  provided 
an  extensive  data  base  of  total  density  measurements  of  the  lower 
thermosphere. 

The  original  accelerometers  were  single  axis  systems  flown 
successfully  on  several  satellites  such  as  the  AE-C,  -D  and  -E  satellites. 
This  instrument  was  called  the  Miniature  Electrostatic  Accelerometer 
(MESA)  and  since  it  was  only  a  single  axis  system,  three  mutually 
perpendicular  sensors  wero  required  to  provide  a  triaxial  sensing 
capability.  The  MESA  system  weighed  10  kg  and  required  20  W  of  power 
(Marcos  and  Swifl,  1982).  The  SETA  instrument  was  an  improvement  over 
the  MESA  by  weighing  less  (4.6  kg)  and  requiring  less  power  (11.7  W). 
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2.4J.  Operation 

As  it's  name  suggests,  the  SETA  instrument  has  a  triaxial  sensing 
capability  as  the  MESA  system,  but  is  accomplished  within  one  sensor.  The 
instrument  is  capable  of  measuring  total  mass  density  and  cross-track 
winds  from  acceleration  measurements  (Marcos  and  Forbes,  1985).  The 
instrument  electrostatically  suspends  a  proof  mass  within  a  hollow 
cylinder  along  its  radial  axis.  The  Z  axis  is  along  the  center  of  the  of  the 
cylinder  and  the  X  and  Y  axes  are  perpendicular  to  the  cylinder.  Each  of 
the  axes  have  a  pair  of  electrodes  that  act  to  constrain  the  proof  mass  (see 
figure  2.2a,b). 

Accelerations  along  each  of  these  axes  are  measured  by  sensing 
motion  in  the  proof  mass  and  applying  an  electrostatic  force  on  the 
appropriate  electrode  that  acts  to  restore  the  proof  mass  to  its  initial  position 
(Marcos  and  Swift,  1982).  In  this  way  the  dc  voltage  required  to  restore  the 
proof  mass  is  proportional  to  the  acceleration.  The  subsequent  output  of  the 
sensor  is  a  digital  pulse  rate  in  proportion  to  the  input  acceleration  (Marcos 
and  Swift,  1982). 

The  accelerations  that  must  be  measured  are  very  small,  so  the 
instrument  sensitivity  range  is  from  1.2  x  10"3  g  to  1  x  10**  g.  However, 
accelerations  as  low  as  2  x  10**  g  can  be  sensed  when  on  the  most  sensitive 
scale  of  the  unit  (Marcos  and  Swift,  1982).  The  instrument  stores  digital 
pulses  for  a  2.045  s  sample  period  and  provides  the  average  acceleration 
over  this  period  as  the  output 

The  instrument  is  calibrated  against  a  portion  of  the  earth’s  gravity 
field  by  tilting  the  instrument  slightly  from  the  horizontal,  zero  position 


(ifiura  22a  Diagram  of  the  Z»axle  Force  Rebalance  Loop  (Marcoe.  1982). 


figure  2.2b  Diagram  of  the  Y-  (or  X-)  axia  Force  Rebalance  Loop  (Marcoe,  1982). 
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(Marcos  and  Swift,  1982).  The  calibration  is  limited  by  a  means  capable  of 
measuring  very  small  changes  in  angular  displacement;  a  laser 
interferometer  is  used  that  allowed  for  a  measurement  as  low  as  10*?g. 

The  instrument  is  mounted  on  the  satellite  so  that  its  Z  axis 
corresponds  to  the  satellites  in-track  direction,  the  Y  axis  is  situated  along 
the  earth’s  radial  direction,  and  the  X  axis  the  cross-track  direction. 

2.4.2  Data 


Typical  plots  of  the  data  represent  X,Y,  and  Z  axis  accelerations  vs 
time  and  is  illustrated  in  figure  2.3.  The  Z  axis  plot  represents  satellite 
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figure  2.3  Filtered  data  from  accelerometer  (Marcos,  1982). 


drag  and  reaches  a  maximum  (~8|ig)  at  the  perigee  where  the  density  is 
highest.  Since  atmospheric  drag  opposes  the  direction  of  the  satellite, 
accelerations  are  negative  and  therefore,  maximum  densities  correspond  to 
more  negative  accelerations  and  minimum  densities  correspond  to  less 
negative  accelerations  (Marcos  and  Swift,  1982), 

Plotted  accelerations  along  the  Y  axis  represent  forcing  radially 
away  from  the  earth  with  spikes  coming  from  the  pitch-attitude  thrusters 
on  the  satellite.  The  plot  for  the  X  axis  shows  the  effects  of  cross  component 
winds  (corresponding  to  zonal  winds  at  low  and  middle  latitudes  and 
meridional  winds  at  high  latitudes  since  the  satellite  is  not  in  a  perfectly 
polar  orbit) . 

The  raw  data  has  contributions  from  vehicle  dynamics,  ay 
(spacecraft  vibrations)  and  instrument  operating  characteristics,  aB  and 
aN  (bias  and  random  noise)  as  well  as  atmospheric  drag  and  wind 

accelerations  (Marcos  and  Forbes,  1985).  The  2.045  s  averaged  instrument 
output  (  aT)  is  therefore  expressed  as 

aTi  =  aw  +  avi  +  att  +  aNi  (2.9) 

where  the  i  represents  the  three  axes,  X,  Y,  and  Z.  Solving  for  aw,  the 
aerodynamic  component 

=  aTi-(aw+aBi  +  aNi)  (2.10) 

we  get  the  acceleration  due  to  aerodynamic  drag  which  is  related  to  total 
atmospheric  density. 

Unwanted  contributions  from  spacecraft  vibrations,  av,  and 
instrument  noise,  aN,  are  removed  by  numerical  filtering  techniques 
(Marcos  and  Swift,  1982).  The  bias  value  aB  is  different  for  each  orbit  and 
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was  removed  by  fitting  a  line  to  average  bias  values  for  22  days  for  each  of 
the  axes  (Marcos  and  Swift,  1981). 

Once  the  corrected  acceleration  term  is  determined,  atmospheric 
density  and  wind  velocities  are  calculated  by  the  expression 

a“ =  ^c‘pVS  (2J1) 

The  values  i,  aw  and  p  are  the  same  as  in  previous  expressions.  is  the 
satellite  frontal  area,  m  is  the  satellite  mass,  C,  are  the  drag  coefficients 
with  respect  to  the  three  axes  velocities.  V  is  the  atmospheric  mass  velocity 
relative  to  the  satellite  and  is  determined  by 

V  =  -V0+VA+VW  (2.12) 

where  V0  ,VA,  and  Vw  respectively  gives  the  velocities  for  the  satellite’s 
inertial  reference  frame,  atmospheric  rotation,  and  neutral  wind  (Marcos 
and  Swift,  1982). 

The  components  for  the  drag  coefficient  are  calculated  using  the 
results  from  Sentman  [1961].  The  drag  coefficients  for  each  axis  is 
calculated  and  plotted  as  a  function  of  pitch  angle  in  figure  11  of  Marcos 
and  Swift  [1982], 

An  iterative  technique  is  used  to  solve  equation  2.11  since  there  are 
more  unknowns  (p,  V)  than  equations.  An  initial  estimate  of  density  is 
determined  and  used  to  calculate  the  three  velocity  components  (Marcos 
and  Swift,  If  32).  The  density  estimate  is  refined  by  using  a  new  V,  and  in 
turn  until  an  acceptable  error  is  reached  (Marcos  and  Swift,  1982). 

Accuracy  of  the  instrument  is  limited  by  accelerometer  bias,  attitude 
uncertainties,  drag  coefficient  estimates,  instrument  noise  not  filtered  out. 
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and  the  inability  to  separate  wind  and  density  effects  along  the  in-track  axis 
(Marcos  and  Swift,  1982).  However,  attitude  errors  are  comparatively 
small,  and  those  due  to  bias  and  noise  are  ~0.5%  at  low  altitudes  (170  km) 
and  climb  by  a  factor  of  five  at  higher  altitudes  (220  km).  The  drag 
coefficients  and  wind  uncertainty  contribute  the  primary  error  in  density 
and  are  within  ±5%  each.  Accelerometers,  therefore,  provide  the  most 
accurate  means  of  measuring  total  density  at  these  low  altitudes. 

2.4.3  Orbital  characteristics 

The  SETA-2  instrument  was  flown  on  a  low  altitude  satellite  with  an 
apogee  of  260  km  and  perigee  of  170  km.  The  data  set  used  in  this  study 
provided  total  mass  density  (in  g/cm^)  measurements  between  May  and 
July  1982.  This  time  period  was  near  solar  maximum  conditions  for  solar 
cycle  21.  The  set  consisted  of  a  measurement  every  20  seconds  totalling 
186,814  data  points  over  1369  orbits. 


figure  2.4  Ft0.7  and  FI 0.7a  values  during  the  SETA-2  period. 
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Figure  2.4  and  2.5  shows  the  F10.7  and  F10.7a,  Ap  and  Kp  indices  for 
the  time  period  of  the  data  set.  The  indices  show  that  there  were  a  number 
of  geomagnetic  disturbances  during  this  time,  providing  ample  data  for  a 
broad  range  of  magnetic  activity  in  order  to  study  this  aspect  of  VSH  and 
MSIS. 

The  orbit  of  the  satellite  took  if  through  an  almost  constant  local  time 
meridian  (10-22  LT)  throughout  the  whole  data  set  with  apogee  being  in  the 
night  sector  and  perigee  in  the  day  sector.  This  obviously  puts  a  local  time 
bias  on  the  results  of  the  study.  However,  the  amount  of  data  available  and 
time  of  coverage  still  provides  a  good  statistical  distribution  to  make 
meaningful  comparisons  with  VSH  and  MSIS. 


120 


40 


2.5  Dynamic  Explorer  NACS  instrument 

The  Dynamics  Explorer  program  was  developed  to  study  interactions 
between  the  magnetosphere,  ionosphere,  and  the  atmosphere.  The 
program  provided  for  two  coplanar  satellites  with  an  inclination  of  90 
degrees  and  with  different  altitudes.  Satellites  DE-1  and  2  were  launched  in 
1981.  The  high  altitude  satellite  (DE-1)  was  used  for  global  auroral 
imaging,  plasmasphere  and  plasmapause  observations,  and  magnetic 
measurements  (Hoffman  and  Schmerling,  1981).  The  low  altitude  satellite 
(DE-2)  had  a  perigee  of  305  km  and  an  apogee  of  1300  km  and  waB  used  to 
provide  neutral  particle  and  wind  measurements.  Orbital  characteristics 
of  the  DE-2  satellite  were  determined  so  that  perigee  went  from  pole  to  pole 
in  less  than  two  months  since  the  emphasis  was  on  polar  phenomena 
(Hoffman  et  al.,1981). 

The  relevant  satellite  for  my  research  is  the  DE-2  satellite. 
Instruments  from  this  satellite  give  neutral  atmosphere  density  and  wind 
data  which  complement  the  lower  altitude  measurements  from  the  SETA-2 
instrument.  By  comparing  measurements  from  these  two  satellites  it  is 
conceivable  to  evaluate  the  height  dependence  of  the  VSH  model. 

One  of  the  nine  instruments  on  this  satellite  is  the  neutral 
atmosphere  composition  spectrometer  (NACS).  This  instrument  measured 
similar  parameters  of  the  neutral  atmosphere,  but  in  slightly  different 
methods  in  order  to  provide  a  redundancy  in  measurements. 


2.5.1  NACS  instrument 
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Density  is  the  NACS  instrument’s  primary  objective.  Instruments 
similar  to  this  one  were  flown  on  th8  AE-C,  -D,  and  -E  satellites  giving  good 
results.  It  uses  a  quadrupole  mass  spectrometer  to  measure  abundances  of 
0,  N2,  He,  and  Ar  in  the  300-500  km  range. 

The  closed  source  mode  is  employed  allowing  the  particles  to  come  in 
chemical  and  thermodynamic  equilibrium.  By  using  this  method, 
chemically  reactive  species  such  as  atomic  oxygen  are  determined 
indirectly  by  measuring  the  species  in  some  reacted  form  (ie.  molecular 
oxygen). 

Particles  enter  the  antechamber  of  the  sensor  where  they  are  allowed 
to  collide  with  the  walls  before  being  ionized  and  measured.  Thi?  enhances 
the  particle  density  as  they  cool  from  1000  K  to  300  K  (temperature  of  the 
sensor)  and  results  in  improved  sensitivity  in  the  detection  of  low 
concentration  particles  such  as  He  and  Ar  in  these  altitude  regions 
(Carignan  et  al.,1981).  The  particles  are  then  ionized  by  electron 
bombardment  in  the  ion  source  region  of  the  sensor  and  are  focused  onto 
the  aperture  of  the  quadrupole.  The  quadrupole  consists  of  four  long  rods 
with  appropriate  potentials  applied  to  allow  mass  selection  with  a 
resolution  of  1  amu.  The  range  for  detection  is  from  2  to  50  amu.  The  ions 
then  have  individual  mass  to  charge  ratios  and  only  ions  with  the  selected 
ratio  move  down  the  rods  to  bo  detected  by  an  electron  multiplier. 

Mass  numbers  are  selected  into  eight  0.016  s  intervals  and  the 
sequence  is  repeated  every  0.1 28  s.  Combinations  of  mass  selection  can  be 
used  to  allow  for  a  particular  investigation  such  as  the  He  to  Ar  ratio.  Mass 
4  could  be  selected  for  two  intervals  and  mass  40  in  the  other  six.  ?V.e  eight 
interval  mass  selection  means  that  a  sweep  of  the  masses  from  1-25  can  be 
accomplished  in  2  seconds  (Carignan  et  al.,1981). 
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Problems  in  instrument  calibration  resulted  in  no  absolute 
calibration  error  being  determined  (A.  Burns,  private  communication). 
Instrument  measurements  drifted  from  month  to  month  with  probable 
error  exceeding  that  of  the  SETA  instrument  measurements,  and  more 
closely  approaching  model  errors.  Subsequently,  comparisons  between 
model  and  NACS  data  are  difficult  to  interpret. 

2i>  Statistical  Techniques 

One  important  aspect  of  the  validation  study  is  the  statistical  method 
ueed  to  validate  the  model  output.  In  choosing  the  method  a  number  of 
items  were  taken  under  consideration.  Obviously,  we  wanted  to  compare 
VSH  density  output  to  the  satellite  measurements,  but  also  wanted  an 
objective  comparison  of  VSH  with  another  model  such  as  MSIS.  We  also 
had  a  large  data  set  to  make  good  statistical  comparisons. 

The  most  simple  way  of  comparing  two  values  is  by  taking  the 
average  relative  error,  or  the  percent  difference  between  them.  This  is 
expressed  mathematically  by 

e  a  fesafatlPga  x  100.  (213) 

P*»U 

We  define  e  to  be  the  percent  difference  of  the  model  (VSH,  MSIS)  and  total 
mass  density  output,  Ptnodel*  relative  to  the  total  mass  density 
measurements  by  the  SETA  instrument,  Psota* 

A  one  to  one  comparison,  using  equation  2.13,  is  taken  between  the 
model  output  and  the  measurements  for  each  point  in  the  SETA  data  set.  A 
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frequency  distribution  is  produced  by  summing  over  the  percent 
differences.  The  resulting  curve  is  a  probability  distribution. 

The  subsequent  plots  generated  illustrate  the  percent  difference 
between  the  models  (VSH  and  MSIS)  and  the  SETA  measurements  versus 
the  respective  frequency  of  occurrence.  Therefore,  two  frequency 
distribution  curves  are  shown,  one  for  VSH  and  the  other  for  MSIS.  Figure 
2.6  illustrates  the  format  for  the  plots  where  the  horizontal  axis  is  the 
percent  difference,  or  “error",  from  SETA  values,  and  the  vertical  axis  is 
the  frequency  (or  number  of  counts)  at  a  specified  error.  The  0  point  is  the 
“true"  (SETA)  value  and  the  areas  of  negative  (positive)  percent  difference  is 
where  the  model  underestimates  (overestimates)  the  SETA  density. 

The  natural  way  of  fitting  a  frequency  distribution  is  by  using  a 
Gaussian  distribution.  The  Gaussian  distribution  is  an  approximation  to 
the  Binomial  distribution  for  the  special  case  where  the  number  of  data 
points  is  very  largo  and  the  probability  of  success  is  very  good  (Bevington, 
1369).  This  certainly  describes  our  case  since  thore  are  almost  200,000  data 
points  and  the  percent  differences,  of  the  model  output,  from  SETA  data  are 
relatively  small. 

Tiie  equation  for  a  Gaussian  distribution  is 


Pq(x,H,o) 


(2.14) 


Figure  2.7  shows  the  symmetric  Gaussian  “bell*  curve  about  a  mean,  p, 
with  standard  deviation,  o.  The  width  of  the  curve,  T,  is  represented  by  the 
full-width  at  half-maximum.  The  tangent  to  the  curve  at  V  drawn  to 
intersect  the  horizontal  axis  locates  the  value  for  two  standard  deviations 
(Bevington,  1369). 
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The  “weighted”  mean  of  the  Gaussian  is  calculated  by 


2>.f< 

\l  =  -  (2.15) 

i-l 

wher8  the  summation  is  over  the  whole  data  set  n,  from  -100  to  +100  percent 
difference  x,  and  the  frequency  f  at  each  xj  (Bevington,  1969).  The 

denominator  is  the  area  bounded  by  the  distribution  curve  and  is  therefore  a 
normalization  factor  to  make  the  sum  of  the  probabilities  equal  to  1.  The 
next  parameter  needed  to  calculate  the  Gaussian  is  the  standard  deviation 
(Bevington,  1969),  which  is  the  square  root  of  the  variance,  o^: 

0!  =  -l=i— ;; - (l1  (2.16) 

I'. 

i-l 


Where  all  values  are  the  same  as  in  equation  2.14. 

The  standard  deviation  is  a  measure  of  the  dispersion  about  the 
mean,  and  is  also  proportional  to  the  uncertainty  in  determining  the  mean 
of  the  distribution  (Bevington,  1969).  Therefore,  the  standard  deviation  is 
the  primary  statistical  value  used  to  make  comparisons  of  VSH  with  MSIS 

One  term  that  is  used  a  number  of  times  to  describe  the  appearance 
of  a  frequency  distribution  curves  is  skewness.  As  defined  by  Press  et  ai. 
[1986],  “skewness  characterizes  the  degree  of  asymmetry  of  a  distribution 
around  its  mean.  It  is  a  pure  number  that  characterizes  only  the  shape  of 
the  distribution.”  Figure  2.8  (Press  et  al.,  1985)  illustrates  positive  and 
negative  skewness,  where  positive  skewness  indicates  a  distribution  with  a 
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tail  that  extends  towards  positive  values  of  x  and  a  negative  skewness 
indicates  a  tail  that  extends  towards  negative  values  of  x.  The  effect  of 
skewness  on  the  plots  in  chapter  3  is  to  shift  the  mean  value  and  sometimes 
increase  the  dispersion  about  the  mean. 


figure  2.8  Graphic  Illustration  of  skewness. 


Problems  were  noted  in  using  a  Gaussian  to  fit  some  of  the  curves. 
In  some  cases  the  VSH  frequency  distribution  was  bifurcated  or  spread  out 
across  large  errors.  The  bifurcation  was  noted  in  an  earlier  version  of  the 
VSH  model  and  the  dispersion  was  noted  when  insufficient  data  existed  to 
compare  with  the  models.  In  these  cases  the  frequency  distribution  was 
distinctly  non-Gaussian  and  the  standard  deviations  calculated  were  not 
good  indicators  to  make  comparisons  between  VSH  and  MSIS. 

However,  this  method  has  been  used  extensively  in  previous  studies 
to  make  statistical  comparisons  (Marcos  et  al.,1989;  Marcos,  1991)  and  it 
was  necessary  to  continue  its  use  in  this  study  to  show  an  historical 
progression  in  comparisons  with  other  models  in  these  studies. 
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In  future  validation  studies  a  more  applicable  statistical  method  is 
recommended.  A  modification  to  the  current  method  would  be  to  weight 
each  of  the  percent  difference  values  based  on  their  difference  from  the 
maximum  frequency  in  the  mean  value  calculation.  A  value  close  to  the 
maximum  value  would  have  a  higher  weight  than  one  further  away.  This 
would  make  the  Gaussian  mean  match  more  closely  with  the  parent 
frequency  distribution  maximum  frequency  value. 
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CHAPTERS 


VALIDATION  STUDY  RESULTS 


3JL  Introduction 

The  previous  two  chapters  provide  the  necessary  background  for  the 
results  of  the  validation  study  included  in  this  chapter.  This  chapter 
includes  a  number  of  plots  that  are  based  primarily  on  SETA-2 
measurements  and  to  a  lesser  extent  on  NACS  measurements  (from  the 
DE-2  satellite)  and  show  certain  dependences  of  the  VSH  model.  NACS 
measurements  were  included  in  an  attempt  to  evaluate  the  altitude 
dependence  of  the  model.  Section  3.2.1  includes  the  results  from  the 
systematic  process  used  to  discover  a  significant  shortcoming  in  the  VSH 
model  and  an  explanation  of  the  subsequent  changes  to  the  model.  Model  to 
model  comparisons  are  then  made  between  the  updated  version  of  VSH  and 
the.  MSIS  empirical  model  in  section  3  2.2. 

&2  Comparison  of  VSH  and  MSIS  to  SETA-2 

A  general  comparison  was  first  made  between  the  models  and  the 
whole  SETA-2  data  set  (with  no  data  binning)  to  get  an  overall  impression  of 
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the  VSH  model’s  performance.  Figure  3.1  illustrates  this  overall  plot.  The 
horizontal  axis  is  the  average  relative  error,  or  percent  deviation  from 
SETA-2  measurements  and  the  vertical  axis  is  the  frequency  or  number  of 
counts.  In  all  of  the  following  plots  the  smooth  dark  lines  are  the  Gaussian 
fits  to  each  of  the  model’s  frequency  distribution  and  the  model  standard 
deviation  and  mean  values  are  indicated  in  the  upper  right  hand  comer. 

It  was  noticed  during  this  study  that  in  a  number  of  cases  the 
Gaussian  fits  do  not  fully  represent  the  frequency  distributions.  Although 
the  Gaussian  fits  are  good  at  showing  model  trends  and  relative 
comparisons  between  the  models,  in  these  cases  the  distributions  are 
simply  non-Gaussian.  However,  the  method  for  calculating  the  standard 
deviation  and  mean  have  been  used  in  past  studies  by  Marcos  [1989]  to  make 
comparisons  between  the  Jacchia  and  MSIS  models.  Therefore,  in  order  to 
make  meaningful  comparisons  between  the  current  VSH  model  and  the 
other  models  in  past  validation  studies  I  have  continued  using  this 
statistical  method. 

3J2L1  Discovery  of  F1Q.7  dependence 

Figure  3.1  is  for  an  early  version  of  VSH  and  was  instrumental  in 
giving  us  our  first  indication  that  a  significant  shortcoming  existed  in  the 
model.  The  VSH  curve  is  distinctly  non-Gaussian,  with  a  prominent  “leg” 
extending  from  -16%  well  into  high  positive  values  (indicative  of  a  positively 
skewed  distribution)  which  indicates  the  model  is  overestimating  density 
under  certain  conditions.  Previous  validation  studies  showed  no  indication 
of  this  skewness  because  the  comparisons  were  made  with  data  over 
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relatively  short  time  periods,  whereas  the  data  set  used  in  this  study  covers 
70  days. 

In  order  to  systematically  investigate  the  performance 
characteristics  of  the  model  and  determine  the  conditions  for  which  the 
VSH  model  is  overestimating  density,  the  SETA-2  data  was  binned  into  low 
(0  <,  Kp  <  3-),  medium  (3-  <,  Kp  <  4+)  and  high  (4+  <.  Kp  <  9)  Kp  indices. 

Figures  3.2  represent  deviations  from  SETA-2  measurements  for  the 
three  magnetic  activity  levels.  The  low  Kp  case  shows  VSH  being  skewed 
towards  positive  values  and  in  the  medium  Kp  case  even  more  skewed  with 
a  larger  area  of  overestimation.  The  VSH  curve  in  the  high  Kp  case  again 
shows  a  wide  area  of  overestimation,  but  now  the  curve  begins  to  show  a 
distinct  bifurcation,  or  secondary  maximum. 

Next,  to  check  for  latitude  dependence,  the  data  was  binned  into  six 
30  degree  latitude  intervals  in  addition  to  the  Kp  bins.  The  top  three  plots  on 
the  page  in  each  of  figures  3.S-3.5  are  for  the  northern  hemisphere  and  the 
bottom  three  plots  are  for  the  southern  hemisphere.  From  loft  to  right  the 
latitude  intervals  indicate  high-latitude  K60°-90°),  middle-latitude  ±{3G°- 
60°),  and  equatorial  ±(0°-30°)  regions. 

Figure  3.3  illustrates  the  latitudinal  dependent  low  Kp  case.  VSH 
continues  to  show  the  positive  skewness,  but  now  the  data  is  becoming  even 
more  bifurcated,  with  small  secondary  and  even  tertiary  maximums 
evident.  In  all  cases,  however,  the  primary  peak  is  very  close  to  the  SETA-2 
measurements  and  shows  little  variation  over  latitude.  The  medium  Kp 
case  in  figure  3.4  shows  similar  results.  The  high  Kp  case  in  figure  3.5 
does  not  illustrate  a  good  statistical  distribution  because  of  the  small 
number  of  occurrences  of  Kp  values  greater  than  4.  But  what  is  most 
evident  in  these  plots  is  how  the  curve  is  being  “stretched"  towards  higher 
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values.  Obviously  the  cause  of  this  bifurcation  is  not  latitude  or  Kp 
dependent  since  it  occurs  through  all  latitudes  and  Kp  values. 

The  root  of  the  problem  was  discovered  when  the  data  was  divided 
even  further  into  F10.7  bins.  Three  F10.7  intervals  were  chosen;  F10.7  < 
160, 160  <i  F10.7  ^  190,  and  F10.7  >190.  Figure  3.6  is  most  indicative  of  the 
results  from  this  binning.  Each  page  represents  one  of  the  latitude 
intervals  with  the  top  row  showing  low  Kp  values,  the  second  row  medium 
Kp  values,  and  the  bottom  row  high  Kp  values.  The  three  columns  on  each 
page  coincide  with  the  three  FI 0.7  bins  mentioned  above. 

The  other  latitude  and  Kp  cases  for  the  FI 0.7  binning  are  located  in 
the  appendix,  a  representative  case  in  included  here  for  evaluation.  The 
series  of  three  FI 0.7  plots  in  figure  3.6  indicates  the  cause  of  the  bifurcation 
in  the  latitude  plot  in  figure  3.4  and  is  for  the  case  of  the  North  equatorial 
region  during  medium  Kp  values.  The  first  F10.7  plot  corresponds  directly 
with  the  primary  peak  of  the  North  equatorial  latitude  plot  of  figure  3.4. 
The  second  FI  0.7  plot  is  an  intermediate  one  that  shows  the  VSH  densities 
beginning  to  drift  towards  higher  values,  while  the  third  FI  0.7  plot  directly 
coincides  with  the  secondary  peak  in  the  latitude  plot.  Therefore,  VSH 
specified  densities  well  during  low  FI  0.7  values  (<  160),  but  began 
overcompensating  as  the  FI 0.7  value  increased  until  a  local  maximum 
occurred  during  high  FI 0.7  values  (>  190). 

This  can  be  further  explained  by  figure  3.7  which  shows  a  plot  of  the 
F10.7  values  versus  the  Kp  values  during  the  SETA-2  period.  The 
horizontal  axis  represents  the  three  F10.7  intervals  used  in  this  study  and 
the  vertical  axis  represents  the  Kp  values  with  horizontal  lines  drawn  to 
represent  the  three  Kp  bins.  This  type  of  plot  illustrates  the  correlations 
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between  the  two  indices  and  helps  explain  that  the  data  was  being 
bifurcated  partly  because  of  the  choice  of  Kp  and  FI  0.7  bins. 

The  figure  illustrates  that  during  low  Kp  values  most  of  the 
occurrences  of  F10.7  values  were  concentrated  in  a  small  interval  (100-160) 
on  the  low  side.  In  the  high  Kp  case,  occurrences  of  FI  0.7  were  mostly  on 
the  high  side,  but  were  spread  out  over  a  larger  interval  (190-280).  The 
narrow  interval  in  the  low  Kp  ca,  ,e  tended  to  concentrate  the  VSH  densities 
into  the  primary  peak  whereas  the  larger  F10.7  interval  tended  to  “spread” 
the  densities  on  the  high  side.  The  sum  of  these  effects  was  to  create  the 
bifurcation. 

The  physical  reason  for  this  problem  was  also  due  to  VSH  using  only 
F10.7  as  the  input  for  solar  flux.  But  F10.7  is  primarily  an  indicator  of  the 
short  term  variability  of  the  sun  which  affects  the  upper  thermosphere 
more  so  than  the  lower  thermosphere  (see  section  1.2. 1.1  for  further 
discussion  of  this  matter).  Therefore  the  SETA-2  measurements  occurred 
at  low  altitudes  this  FI 0.7  dependence  showed  up  indicating  a  need  for  an 
additional  parameterization.  The  F10.7a  index  was  included  in  a  height 
dependent  solar  flux  parameterization  (section  1 .2.1.1)  since  it  better 
indicates  the  state  of  the  lower  thermosphere,  The  current  version  of  VSH 
(with  this  parameterization)  shows  marked  improvement  in  density 
specification  and  shows  no  bifurcation  of  data.  Results  of  the  current 
version  are  included  in  section  S.2.2. 


&2.U  Height  variation 

Another  problem  noted  during  this  study  has  to  do  with  the  NCAR- 
TIGCM  and  its  ability  in  specifying  the  mean  density  values.  The  TIGCM 


Correlation  between  the  Kp  bins  and  the  F10.7  bins  used 
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is  very  good  at  computing  variations  of  output  fields  about  a  global  mean. 
However,  as  mentioned  previously,  in  the  case  of  neutral  density  output,  the 
TIGCM  global  mean  is  different  from  that  of  the  MSIS  model,  which  is 
derived  from  climatologically  averaged  satellite  data  (Killeen,  et  al.,  1991). 
The  density  mean  of  the  TIGCM  has  to  be  adjusted  to  be  more  in  line  with 
MSIS. 

It  was  initially  thought  that  the  density  profile  of  the  for  the  TIGCM 
and  MSIS  looked  like  figure  3.8a,  where  the  scale  heights  (or  slopes)  were 
equal  but  offset  by  a  factor.  Based  on  this,  one  height  level  was  chosen  in 
making  the  correction  factor  between  the  two  scale  heights  (a  to  b  in  figure 
3.8a).  But  later  it  was  found  that  the  scale  heights  were  different  as  in 
figure  3.8b.  A  correction  was  made  in  going  from  a  to  b  and  then  also 
applied  to  a  lower  level,  following  the  assumed  scale  height,  the  density  was 
still  off  by  the  difference  between  c  and  d.  Therefore,  a  single  height  level 
couldn’t  provide  a  proper  correction  factor.  A  second  height  was  used  to  get 
an  improved  adjustment  as  in  figure  3.8c  and  correction  factors  were 
determined  by  interpolating  between  the  two  heights. 

This  adjustment  is  important  in  mapping  the  proper  height  to  a 
pressure  level  (which  is  the  vertical  coordinate  system  used  in  the  TIGCM) 
ior  use  by  the  VSH  model.  This  also  results  in  enhancing  the  performance 
of  the  mass  density  output  from  VSH. 

3.2.2  Model  to  data  comparisons  using  current  version  of  VSH 

The  results  included  in  this  section  come  from  the  August  1991 
version  of  VSH.  I  used  the  same  systematic  approach  with  the  current 
version  as  with  the  older  version,  showing  an  overall  dependence  of  the 
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figure  3.8a  Scale  height  comparison  (expected) 


Density 


figure  3.6P  Scale  height  comparison  (expected) 


Density 

fig.  3  8c  Scale  Heigh  Comparison  (actual) 
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model  to  the  SETA-2  data  and  then  binning  the  data  into  Kp  and  latitude 
intervals.  In  addition  to  this  I  have  also  included  MSIS  curves  to  make 
inter-model  comparisons  and  also  identify  and  attempt  to  explain  the 
physical  processes  that  cause  model  variability.  Since  the  standard 
deviation  value  is  most  indicative  of  how  well  a  model  specifies  density,  I 
will  primarily  use  this  value  as  a  “yardstick”  to  compare  VSH  and  MSIS 
results. 

The  last  subsection  includes  model  to  model  comparisons  of  standard 
deviation  results  based  on  binning  by  various  parameters  including  Kp, 
latitude,  altitude,  and  FI  0.7.  These  results  indicate  that  VSH  does  as  well 
as  and  in  many  cases  significantly  better  than  MSIS  in  the  specification  of 
density, 

3&2L1  Overall  dependent 

The  whole  data  set  was  compared  with  both  VSH  and  MSIS  models 
with  no  binning  of  data.  The  results  of  this  are  illustrated  in  figure  3.9.  As 
can  be  seen  the  VSH  curve  is  markedly  improved  over  tiro  older  version  of 
VSH  (see  figure  3.1  for  comparison).  The  frequency  distribution  for  VSH  is 
more  Gaussian-likc  with  little  dispersion  about  the  mean  and  no  areas 
where  the  density  is  strongly  over-  or  underestimated.  In  comparison,  the 
MSIS  curve  has  more  dispersion  about  the  mean,  creating  a  larger 
standard  deviation  than  that  of  VSH  (14.8%  compared  to  12.2%  for  VSH). 
This  standard  deviation  for  MSIS  corresponds  well  with  previous  studies 
(Marcos, 1989)  and  represents  the  limit  of  performance  for  empirical 
models.  In  comparing  the  mean  values,  MSIS  is  much  closer  to  the  SETA- 
2  measurements,  whereas  the  VSH  mean  indicates  underestimation  by 
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figure  3.9  OwaU  plot  of  now  VSH  deviations  from  SETA  values. 
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7.0%.  However,  for  the  end  user  of  the  VSH  model  a  simple  correction 
factor  can  be  included  to  shift  the  mean  in  the  desired  direction. 

Little  physical  information  can  be  derived  from  this  plot,  only  to  note 
that  MSIS  is  slightly  skewed  to  the  right  with  a  small  “tail”  where  a 
number  of  density  values  are  being  overestimated  by  30-40%.  However,  no 
reasoning  for  this  can  be  given  based  on  this  plot.  More  information  may  be 
concluded  by  dividing  the  data  under  various  conditions. 

3. 2.2.2  Kp  dependence 

Using  the  same  procedure  as  section  3.2.1  the  data  is  binned  into 
three  Kp  intervals  in  order  to  determine  the  performance  characteristics  of 
the  model  based  on  magnetic  activity.  Figure  3.10  illustrates  low,  medium 
and  high  Kp  conditions  for  VSH  and  MSIS. 

The  models  represent  density  best  during  low  Kp  cases.  The  VSH 
curve  has  little  dispersion  and  does  not  appear  to  be  skewed  in  either 
direction.  MSIS,  on  the  other  hand,  has  a  larger  standard  deviation  than 
VSH  but  a  better  handle  on  the  mean  value.  As  in  figure  3.9  MSIS  is 
positively  skewed  indicating  some  conditions  where  densities  are  being 
overestimated. 

Under  the  medium  Kp  case  it  can  be  seen  that  both  models  begin  to 
disperse  more  and  the  mean  values  are  shifted.  In  the  high  Kp  case  (note 
change  in  vertical  scale)  the  models  have  a  more  difficult  time  representing 
density.  The  standard  deviations  for  each  increase  to  values  >16%.  Both 
models  are  also  strongly  skewed  positively,  more  so  in  the  MSIS  curve, 
which  accounts  for  the  larger  standard  deviation. 
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The  MSIS  model  does  poorer  in  the  high  Kp  case  than  in  low  Kp  case,  due  to 
the  spatial  resolution  of  the  model.  As  described  in  section  2.2,  MSIS 
retains  ~5  harmonics  which  in  turn  limits  the  effective  spatial  resolution  to 
20  degrees.  The  VSH  model  also  has  poorer  statistics  during  high  Kp 
cases,  but  tends  to  do  slightly  better  (a  standard  deviation  of  16.2% 
compared  to  19.7%  for  MSIS)  possibly  due  to  the  better  resolution  of  the 
NCAR-TIGCM,  which  has  an  effective  spatial  resolution  of  10  degrees.  The 
limited  spatial  resolution  of  the  models  makes  it  difficult  to  resolve  the  fine 
density  structure  that  occurs  during  times  of  high  Kp  such  as  during 
geomagnetic  storms.  During  these  storms  the  density  can  vary  greatly 
from  the  high  density  area  in  the  “throat”  of  the  ion  convection  pattern  to 
the  relatively  low  density  surrounding  areas. 

The  mean  values  for  both  models  drift  from  negative  to  more  positive 
values  from  low  to  high  Kp  values.  This  is  most  probably  due  to  the  models 
overcompensating  as  the  Kp  levels  increase.  This  pattern  of  the  mean 
drifting  is  also  obvious  in  the  latitude  plots  shown  in  the  next  subsection. 

3uL2L3  Latitude  dependence 

The  accuracy  of  thermospheric  models  is  limited  by  their  ability  to 
resolve  atmospheric  occurrences  that  happen  on  very  small  scales.  As 
described  in  the  first  chapter  the  ion  convection  pattern  in  the  high  latitude 
regions  causes  intricate  circulation  patterns,  upwelling  and  enhanced 
localized  heating  that  results  in  small  scale  variability  in  the  density 
structure.  During  geomagnetic  storms  when  the  Kp  index  increases,  the 
variability  in  density  becomes  even  more  complicated  since  the  circulation 
pattern  is  even  further  enhanced  and  distorted. 
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The  next  step  in  evaluating  model  performance  was  to  put  the  data  in 
six  30  degree  latitude  bins  while  keeping  the  Kp  bins.  The  plots  in  figures 
3.11-3.13  are  set  up  as  described  in  section  3.2.1. 

The  VSH  and  MSIS  curves  exhibit  little  dispersion  in  the  low  latitude 
regions  during  low  Kp  values,  with  VSH  having  slightly  better  standard 
deviations,  but  both  curves  become  more  dispersed  in  the  middle  and  high 
latitudes.  Also  noticed,  is  that  both  models  show  positive  skewness  with 
increasing  latitude,  possibly  due  to  the  models  beginning  to  compensate  for 
the  more  structured  auroral  regions  but  tend  towards  overestimation. 

In  the  medium  Kp  case  both  models*  standard  deviations  increase  by 
-1-2%  over  the  low  Kp  case  due  to  the  enhanced  density  structure  that 
begins  as  magnetic  activity  increases. 

As  seen  in  the  low  Kp  case  the  models  do  well  in  the  low  latitudes  but 
their  standard  deviations  increase  by  a  factor  of  two  in  the  high  latitude 
case.  The  models  begin  to  show  this  increase  in  the  mid-latitudes  due  to 
auroral  effects  becoming  evident  even  at  mid-latitudes.  As  mentioned 
before,  the  resolving  capability  of  the  models  limit  their  ability  to  specify  the 
fine  density  structure  that  shows  up  during  those  effects. 

It  is  also  evident  that  there  is  some  pattern  to  the  way  the  mean  value 
drifts.  This  is  evident  in  each  of  the  cases  to  some  degree.  It  is  most 
obvious  in  the  medium  Kp  case  of  the  MSIS  curve.  The  moan  value 
increases  from  negative  to  more  positive  values  as  the  latitude  goes  from 
the  southern  to  the  northern  hemisphere.  This  occurs  to  a  lesser  extent  in 
the  VSH  curves  also, 

There  really  is  not  enough  information  to  definitively  explain  the 
occurrence  of  moan  drift.  It  could  be  due  to  a  seasonal  dependence  of  the 
model,  with  underestimation  occurring  in  the  southern  hemisphere  winter 
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and  overestimation  occurring  in  the  northern  hemisphere  summer.  Or,  it 
could  simply  be  due  to  hemispheric  differences.  Pole  to  pole  differences  in 
the  turbopause  height  would  allow  for  differences  in  total  density 
measurements.  Also,  since  the  circulation  pattern  is  slightly  different  in 
the  southern  high  latitudes  than  in  the  northern  high  latitudes,  resulting 
density  measurements  could  also  be  different. 

Another  feature  evident  in  these  plots  is  that  VSH  tends  to  do  better  in 
the  northern  hemisphere  than  in  the  southern  hemisphere,  whereas  MSIS 
shows  no  hemispheric  bias.  This  bias  of  VSH  is  possibly  due  to  the 
northern  hemisphere  bias  exhibited  in  the  NCAR-TIGCM. 

In  the  high  Kp  case,  degradation  in  the  model  curves  is  quite  evident, 
especially  in  the  high  latitude  region.  VSH  is  strongly  skewed  in  the 
positive  direction  in  the  southern  hemisphere,  but  MSIS  is  even  more  so  in 
all  cases  with  a  secondary  maximum  beginning  to  show  up  (see  south  mid¬ 
latitude  and  equatorial  region  plots).  This  secondary  maximum  is  most 
possibly  a  high  solar  flux  dependence  by  MSIS.  The  plot  in  figure  3.14 
illustrates  the  same  FI 0.7  bins  as  in  section  3.2.1  and  shows  that  the  area  of 
overestiraation  in  the  MSIS  curve  does  not  show  up  in  the  low  F10.7  plot  but 
does  appear  to  occur  in  the  high  FI 0.7  plot. 

The  next  series  of  plots  show  the  results  from  comparing  the  models 
to  measurements  from  the  NACS  instrument  on  the  DE-2  satellite  during 
the  same  time  period  as  the  SETA-2  data  set.  This  satellite  flew  at  an 
altitude  range  of  300  to  over  1000  kra  with  composition  measurements  were 
taken  only  up  to  -500  km.  It  was  hoped  that  using  high  altitude  DE-2  data 
in  conjunction  with  the  low  altitude  data  from  the  SETA-2  instrument 
would  help  show  altitude  dependences  of  the  models.  However,  upon 
comparing  the  NACS  plots,  figures  3.15-3.17,  to  the  SETA-2  plots  it  became 
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quite  obvious  that  the  absolute  instrument  accuracy  approaches  model 
accuracy  due  to  the  calibration  problem  explained  in  section  2.5.  Also,  the 
amount  of  NACS  data  during  this  time  period  was  about  a  factor  of  ten  less 
than  there  was  SETA-2  data,  and  the  resulting  statistical  distributions  (note 
vertical  scale)  are  not  quite  as  good.  This  small  amount  of  data  and  the 
instrument  accuracy  make  model  to  data  comparisons  quite  difficult  to 
make. 

Absolute  comparisons  cannot  be  made  with  the  SETA-2  plots  for  the 
above  reasons  but  comparisons  can  be  made  and  trends  can  be  detected 
among  the  two  models  within  the  data  set.  As  in  the  previous  section 
figures  3.15-3.17  represent  latitude  and  Kp  dependences.  There  was  more 
northern  hemisphere  data  than  southern  hemisphere  data  available 
during  this  time  since  perigee  was  in  the  northern  high  latitudes,  so 
Gouthem  latitude  comparisons  are  difficult  to  make  and  are  not  included  in 
the  following  short  discussion. 

In  the  low  Kp  case  MSIS  shows  better  standard  deviations  than  VSH 
over  all  latitudes,  but  the  mean  values  are  fairly  close  to  each  other.  The 
VSH  mean  values  show  a  definite  trend  towards  higher  errors  with 
increasing  latitude  for  all  three  Kp  cases.  This  could  show  a  height 
dependence  of  the  model  or  it  could  be  variability  in  the  instrument 
accuracies  over  altitude  ranges.  Little  else  can  be  said  about  the  medium 
and  high  Kp  cases  except  that  the  standard  deviations  of  both  models  are 
close  through  all  latitudes. 
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&SL2.4  Standard  deviation  comparisons 

To  make  easily  observable  comparisons  between  VSH  and  MSIS  I 
binned  the  SETA-2  data  into  four  parameters  separately.  The  standard 
deviations  of  VSH  and  MSIS  were  calculated  for  each  case  and  compared 
using  a  bar  graph. 

The  parameters  include  Kp,  latitude,  FI 0.7,  and  altitude.  The  first 
three  parameters  were  binned  using  the  same  intervals  as  before  and  the 
altitude  parameter  was  binned  into  three  intervals  over  the  SETA-2  altitude 
range  (160-200, 200-230, 230-260). 

Figure  3.18  shows  the  results  from  binning  the  data  by  low,  middle 
and  high  Kp  values.  In  each  case  VSH  has  a  lower  standard  deviation  than 
MSIS  by  ~2  percentage  points. 

The  next  figure,  3.19,  indicates  the  variability  of  the  standard 
deviation  over  latitude.  VSH  does  very  well  in  the  northern  hemisphere  as 
compared  to  MSIS,  but  the  VSH-MSIS  difference  is  much  closer  in  the 
southern  hemisphere  with  higher  standard  deviations  for  each 
corresponding  latitude  interval  in  the  southern  hemisphere.  This  is  most 
probably  due  to  the  hemispheric  bias  of  the  NCAR-TIGCM  as  mentioned 
earlier  in  this  chapter. 

The  F10.7  standard  deviations  of  the  models  in  the  next  figure  are 
quite  close  in  the  medium  and  high  FI  0.7  cases  but  VSH  does  ~3  percentage 
points  better  than  MSIS  in  the  low  FI  0.7  case. 

The  final  graph  shows  the  altitude  dependence  of  the  models.  Again, 
VSH  does  significantly  better  (~4%)  in  the  lower  altitude  regions,  but  the 
two  models  are  essentially  equal  in  the  high  altitude  interval. 
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fiflure  3.20  Same  as  figure  3.18  except  for  F10.7  dependence. 
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figure  3.21  Seme  as  figure  3,18  except  for  attitude  dependence. 
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These  results  were  determined  using  the  same  statistical  technique 
described  in  chapter  1  and  are  consistent  with  the  results  from  the 
validation  study  in  the  previous  subsection. 

3*3  Conclusions 

This  is  one  of  the  first  extensive  validation  studies  of  the  VSH  model. 
Neutral  density  output  from  the  model  was  compared  to  accelerometer 
measurements  of  the  lower  thermosphere.  The  SETA-2  data  set  was  quite 
extensive,  providing  70  days  worth  of  data  during  solstice  conditions  in  1982 
which  was  a  time  close  to  solar  maximum.  The  SETA-2  measurements 
were  shown  to  be  quite  helpful  in  validating  the  VSH  model  in  a  region  of 
the  atmosphere  that  is  quite  data  sparse. 

More  work  can  be  done  using  this  data  set,  such  as  evaluating  any 
local  time  dependences  of  the  model.  After  this,  the  validation  study  can 
continue  by  using  measurements  from  other  instruments,  such  as  the 
SETA-1  instrument  that  took  lower  thermospheric  wind  and  total  density 
measurements  during  equinox  conditions  of  1979. 

During  this  study  we  discovered  a  significant  shortcoming  in  the 
VSH  model  based  on  a  solar  flux  dependence.  The  problem  was  corrected 
using  a  height  dependent  PI 0.7  parameterization  that  includes  the  use  of 
the  FI 0.7  average  value.  The  subsequent  results  show  marked 
improvements  in  density  specification. 

The  data  was  binned  using  a  variety  of  parameters  to  test  the 
characteristics  of  the  VSH  and  MSIS  models.  The  results  show  that  VSH 
exhibit  very  low  standard  deviations  (<  10%)  during  low  magnetic  activity 
and  virtually  all  latitude  regions.  Only  minor  increases  of  ~5%  for  each 
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latitude  region  are  shown  as  the  magnetic  activity  increases  to  higher  Kp 
values. 

In  the  overall  comparison,  with  no  data  bins,  we  found  that  VSH 
exhibited  a  standard  deviation  of  12.2%  or  ~3%  lower  than  the  MSIS 
empirical  model.  Also,  during  tho  model  comparisons  based  solely  on  Kp, 
latitude,  FI  0.7  and  altitude,  VSH  standard  deviations  values  were 
noticeably  less  in  all  cases  than  MSIS  values.  In  particular,  the  northern 
high  latitude  and  low  altitude  regions  are  -3-4%  lower. 
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CHAPTER  4 


EFFECTS  OF  ORBITAL  IN-TRACK  WINDS  ON  SATELLITE  DRAG 

4JL  Introduction 

Neutral  winds  are  important  in  the  determination  of  the  forces 
experienced  by  satellites.  Such  winds  in  the  polar  cap  region  of  the  earth's 
thermosphere  can  reach  speeds  greater  than  1000  m/s  and  can  change  the 
drag  or  lift  forces  felt  by  a  satellite  by  a  largo  amount  duo  to  the  uncertainty 
in  the  wind  velocity.  Such  changes  can  subsequently  affect  the  orbit  of  the 
satellite.  Also,  instruments*  such  as  SETA,  measuring  density  are 
dependent  on  wind  information  for  being  able  to  accurately  specify  the 
density  (Killeen  et  ah,  1983).  Obviously,  it  is  important  to  be  able  to 
accurately  measure  winds  and  know  the  effects  that  such  winds  have  on 
space  structures. 

In  this  chapter  I  have  included  the  results  from  my  research  that 
attempted  to  quantify  the  effects  of  meridional  winds  on  a  polar-orbiting 
satellite.  My  goal  was  to  conduct  this  study  over  a  large  number  of  different 
orbits  in  order  to  provide  a  good  statistical  distribution  of  this  effect.  I  also 
wanted  to  calculate  the  integrated  drag  through  a  whole  orbit  and 
determine  the  percentage  of  the  drag  over  different  regions  about  the  globe. 
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However,  due  to  time  constraints,  only  the  drag  effect  at  each  orbital  point 
was  determined.  Rather  than  using  measured  winds,  I  used  VSH  model 
neutral  meridional  winds  to  illustrate  the  research  capabilities  of  the 
model.  Initial  results  indicate  that  a  700  m/s  wind  parallel  to  the  satellite 
velocity  vector  changes  the  satellite  drag  by  23%  during  a  time  of  high 
magnetic  activity. 

42  Procedure 

To  begin  with  I  chose  a  typical  DE-2  orbit  late  in  the  satellite’s 
operational  cycle  in  order  to  get  a  more  circular  orbit.  A  closely  circular 
orbit  would  allow  primarily  latitude  dependent  wind  variations  and  cut  out 
most  of  the  altitude  dependence.  The  altitude  interval  through  the  orbit  was 
300  *  500  km  with  the  perigee  over  the  high  latitude  regions.  The  orbit  was 
also  constrained  to  certain  local  time  meridians  since  we  wanted  the 
satellite  to  pass  parallel  to  the  transpolar  flow,  or  the  middle  of  the 
circulation  pattern.  The  center  of  tho  circulation  pattern  lies  roughly 
parallel  to  the  0-12  local  time  meridian  and  varies  widely  during  enhanced 
magnetic  activity. 

In  order  to  provide  for  a  worse-case  scenario,  a  time  was  chosen 
during  1982  when  the  FI  0.7  and  Ap  indices  were  high  indicating  strong 
magnetic  activity  and  an  enhanced  ion  convection  pattern.  July  14  was 
selected  with  the  following  geophysical  conditions: 

F10.7 . 269.2 

BT0.7a....l72.7 
Ap . 153 
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The  next  step  included  using  a  FORTRAN  routine  that  created  the 
spatial  and  temporal  points  along  the  orbit  given  the  initial  conditions  of 
apogee,  perigee  and  local  time.  Along  with  the  indicated  geophysical 
conditions  these  orbital  points  were  used  by  the  VSH  computer  model  to 
provide  the  total  mass  density  and  the  neutral  meridional  winds  required 
for  satellite  drag  calculations.  The  method  used  for  this  calculation  is 
described  in  section  1 .3. 

The  resulting  output  was  used  in  an  IDL  routine  that  determined  the 
aerodynamic  drag  based  on  the  in-track  (meridional)  winds  and  the  density 
along  the  orbit.  The  output  from  this  routine  included  satellite  drag  due  to 
winds  and  the  satellite  drag  without  winds  being  taken  into  account.  Using 
these  two  values  the  change  in  drag  effect  due  to  winds  was  calculated  for 
each  orbital  point. 

4.3  Results  and  discussion 

The  data  illustrated  in  tins  section  results  from  an  orbit  with  a  02-1 4 
local  time  meridian  with  the  geophysical  conditions  specified  in  the 
previous  section.  After  observing  the  meridional  wi.ids  along  orbits  with 
different  meridians  it  was  determined  that  this  local  time  slice  exhibited 
the  strongest  winds  (742  m/s  at  a  latitude  of  -  85°)  over  tiro  southern  polar 
cap  region. 

The  percent  difference  between  the  aerodynamic  drag  duo  to  winds 
and  the  drag  experienced  without  the  effect  of  winds  was  calculated  and 
plotted  against  the  latitude  along  the  orbit.  The  results  are  shown  in  figure 
4,1 .  The  plot  shows  that  the  largest  absolute  value  of  the  percent  difference 
causes  a  change  in  the  aerodynamic  drag  by  23%  due  to  a  wind  speed  of  742 
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Latitude 

figure  <1  Percent  difference  between  dreg  dye  to  winds  and 
dreg  wiUi  no  winds  plotted  i&titudinally 
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m/s  over  the  southern  pole.  Since  the  percent  difference  is  a  negative  value 
(the  meridional  wind  vector  is  parallel  to  the  satellite  velocity  vector)  this 
actually  means  that  the  drag  decreases  or  lift  increases  by  23%.  However, 
due  to  the  symmetry  of  the  conditions,  if  the  wind  was  directed  against  the 
satellite  velocity  vector,  the  drag  would  increase  by  23%. 

The  plot  coincides  directly  with  the  direction  of  the  meridional  winds 
but  not  with  the  magnitude,  so  a  positive  percent  difference  coincides  with  a 
northerly  directed  wind  and  the  opposite  goes  for  the  negative  values.  It  is 
then  interesting  to  note  that  the  plot  also  identities  the  circulation  pattern 
exhibited  over  the  polar  cap  region.  This  is  shown  by  the  sudden  change  in 
the  sign  of  the  values  as  the  plot  goes  from  positive  to  negative  and  back  on 
either  side  of  the  large  peak  in  the  southern  hemisphere.  A  large  peak  is 
also  noticed  in  the  northern  high  latitudes  which  coincides  with  the 
circulation  pattern  in  this  region. 

Although  the  results  of  this  truncated  study  are  inconclusive,  they 
still  indicate  that  winds  can  play  a  major  part  in  determining  satellite 
drag,  especially  during  strong  geomagnetic  activity.  This  study  also 
showed  that  Urn  VSH  model  is  quite  useful  in  conducting  research  that 
requires  input  of  thermospheric  parameters.  Rather  than  attempting  to 
search  a  large  satellite  database  for  a  particular  orbit  with  the  required 
orbital  characteristics  and  conditions  when  a  certain  instrument  is  turned 
on,  I  was  able  to  acquire  the  information  much  quicker  and  easier  with  the 
VSH  model. 
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CHAPTERS 


DISCUSSION 

A  validation  study  of  the  total  mass  density  output  from  the  Vector 
Spherical  Harmonic  (VSH)  computer  model  has  been  conducted  with  the 
subsequent  results  included  in  this  thesis.  Results  from  a  small  amount  of 
research  conducted  on  the  effects  of  VSH  neutral  in* track  winds  on  satellite 
drag  have  also  been  'provided.  Additional  background  on  the  thermospheric 
region  and  aerodynamic  drag  was  included  at  the  beginning  of  the  thesis.  X 
also  provided  information  on  the  tools  used  in  this  study,  which  include  the 
VSH  and  MSXS  models,  data  from  the  SETA-2  and  NACS  instruments,  and 
finally  Urn  statistical  method  used  to  compare  the  results  in  an  objective 
manner. 

An  extensive  validation  program  now  exists  that  can  to  used  to 
compare  the  VSH  model  to  any  data  set.  The  program  consists  of  a 
FORTRAN  program  that  ht  capable  of  binning  the  data  by  various 
parameters  and  comparing  these  measurements  to  the  corresponding 
results  from  the  VSH  and  MSXS  models.  The  output  from  this  is  then  ran 
through  another  FORTRAN  program  that  calculates  the  mean  value  and 
standard  deviation  of  tho  specified  data  field  and  fits  a  normal  (Gaussian) 
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distribution  to  the  data.  An  IDL  program  allows  the  user  to  specify  the 
particular  bins  and  then  plots  the  results  in  an  appropriate  format. 

The  plots  illustrating  the  magnetic  activity  and  latitude  dependence 
of  the  models  show  that  VSH  provides  the  best  results  (with  standard 
deviations  less  than  10%)  during  times  of  low  Kp  values  and  in  lower 
latitude  regions.  There  is  only  a  factor  of  1/2  increase  during  high  Kp 
values  and  in  higher  latitudes.  The  VSH  model  does  tend  to  show  a 
hemispheric  bias  with  better  results  in  the  northern  latitudes. 

Finally,  in  an  overall  comparison  (through  the  whole  data  set  with  no 
bins)  of  VSH  and  MSIS,  the  standard  deviation  for  VSH  is  ~3  percentage 
points  lower  than  that  of  MSIS.  A  value  of  12.2%  fcr  VSH  is  much  lower 
than  the  standard  that  has  been  exhibited  by  empirical  models  over  the  past 
decade.  This  shows  that  a  mode)  based  or.  truncated  output  of  a  TGCM  is 
capable  of  improving  our  ability  to  specify  thermospheric  density. 

However,  it  must  be  mentioned  that  these  results  are  based  solely  on 
the  SETA-2  measurements.  Although  an  accelerometer  provides  some  of 
our  most  accurate  measurements  ef  total  density  in  the  lower 
thermosphere,  the  information  provided  by  the  SETA-2  measurements  is 
limited  by  the  orbital  <  haracterictics  of  the  satellite  it  flew  on.  There  was  a 
strong  local  time  bias  of  the  data,  all  the  measurements  are  at  an  almost 
constant  local  time  meridian  (10-22  LT),  Also,  perigee  occurred  in  the  day 
sector  and  apogee  in  the  night  sector. 

Another  limitation  to  the  results  includes  the  statistical  method  used 
in  comparing  the  two  models.  Although  the  technique  for  fitting  a 
Gaussian  distribution  to  the  data  is  the  same  as  in  previous  studies,  due  to 
the  extent  of  the  data  and  the  extent  of  the  skewness  exhibited  by  the 
distributions,  many  of  the  model  frequency  distributions  were  distinctly 
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non-Gaussian  and  the  use  of  these  statistics  did  not  fully  represent  the 
distribution.  The  resulting  fits  shifted  the  mean  value  and  showed  more 
dispersion  about  the  mean  than  the  main  body  of  the  distribution. 
Therefore,  more  robust  statistics  are  recommended  in  future  validation 
studies.  Such  a  technique  would  strongly  weight  values  closer  to  the  mean 
than  those  that  are  further  away  and  then  use  a  full  width  at  half 
maximum  method.  This  method  would  apply  only  to  the  maip  body  of  the 
distribution  rather  than  taking  into  account  the  “outliers’*,  or  values  further 
away  from  the  mean. 

Nonetheless,  with  the  SETA-2  conditions  taken  into  account  VSH  still 
provides  very  good  results.  Under  most  conditions  the  VSH  curves  showed 
little  dispersion  about  the  mean  and  exhibited  little  or  no  latitude  dependent 
drift  in  the  mean  value.  There  were  only  a  few  situations  where  the 
Gaussian  fits  were  not  representative.  In  those  cases,  the  VSH  distribution 
curve  was  nevertheless  quite  close  to  the  MS1S  curve. 

The  results  from  the  research  conducted  on  the  effects  of  wind  drag 
on  satellites  are  not  dramatic,  but  do  tend  to  support  previous  studies 
indicating  the  importance  of  neutral  winds  on  satellite  drag.  The  initial 
results  indicated  a  23%  change  in  the  drag  experienced  by  the  satellite  due 
to  a  740  m/s  wind.  Since  winds  in  the  polar  cap  region  have  been  measured 
at  much  higher  speeds  during  strong  magnetic  activity  it  is  obvious  that 
the  drag  felt  by  a  satellite  flying  through  such  an  area  will  be  even  more 
changed  due  to  the  velocity  squared  term  in  the  equation  for  aerodynamic 
drag. 

The  gist  of  this  study  has  to  do  with  oui  ability  to  specify  the 
parameters  that  effect  satellites  in  orbit  about  the  earth.  Both  winds  and 
density  have  a  strong  affect  on  such  satellites  and  our  ability  to  maintain 
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and  use  these  spacecraft  can  be  strongly  impacted  by  the  accuracies  with 
which  we  are  capable  of  specifying  these  parameters.  Further  research 
into  such  computer  models  as  VSH  will  only  increase  our  knowledge  of  the 
thermosphere  and  allow  us  to  better  keep  track  of  the  valuable  resources 
that  orbit  through  this  region. 
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APPENDIX  A 
F10.7  dependence  plots 

This  appendix  contains  the  plots  described  in  section  3.2.1  that 
illustrated  the  PI  0.7  dependence  depicted  by  the  older  version  of  VSH.  The 
plots  on  each  page  represent  one  latitude  region.  The  three  columns 
represent  the  three  FI  0.7  intervals  and  the  rows  represent  the  three  Kp 
intervals.  These  plots  show  conclusively  that  the  VSH  curves  were 
biiurcating  based  on  solar  flux  inputs;  ie.  the  large  positively  skewed  area 
on  the  latitude  plots  occur  during  high  FI  0.7  values  as  indicated  by  the  plots 
in  this  appendix. 
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